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ABSTRACT 


A Navier-Stokes problem solver, developed by L. N. Sankar, 1s modified to provide 
dvnamuc, interactive graphical presentations of predicted flow field solutions about a 
NACA-0012 airfoil section, oscillating in pitch, tn order to demonstrate the capabilities 
of dynamic graphics applications in the study of complex, unsteady flows. Flow field 
solutions in the form of pressure coefficient and stream function contour plots about an 
airfoil experiencing dynamic stall are plotted utilizing an IRIS 3000-series workstation 
and Graphical Animation System (GAS) software, developed by Sterling Software for 
NASA. These full cvcle solutions, in conjunction with dynamic surface pressure dis- 
tribution plots and integrated lift, pitching moment and drag coefficient data, are com- 
pared to existing experimental data in order to provide an indication of the validity of 
the code’s far-field solution. Full procedural documentation is maintained in order to 
provide an efficient analysis tool for use in future oscillating airfoil studies planned by 
the NASA-Ames Fluid Mechanics Laboratory and the Naval Postgraduate School De- 


partment of Aeronautics and Astronautics. 
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1. INTRODUCTION 


It has been recognized for some time [Rel. E] that an airíoil oscillaung in pitcli to 
angles of attack greater than the static stall angle will surpass the traditional stall barier 
e mormal forces whiclr exceed those alttuinable 1 the suatic case. This dv- 
namuc stall phenomenon is attributed to the alt inovement ofa strong shed vortex elong 
и ррег зиасе of the antodl, carrume withit au induced velocity field which radically 
e annal changes eherdwise pressure distributions. In general. an accurate in- 
terpretauou of the dynamic stall mechanisin will siguiftcantly impact a variety cf appii- 
cations, all of which involve dynanue lifting surface motions and unsteady [low 
separation. Originally, dynanuc stall analysis efforts were directed toward helicopter 
aerodvnamics, where sharp increases in oscillatory torsional loading and thus, blade 
stress, can reduce the fatigue Irfe of rotor mechanical components and vortex-induced 
aerodvnamic loading can geneiate adversely phased pitclung moments. resulting in stall 
flutter [Ref. 2]. Much current interest concerns tlie feasibility of exploiting dvnamue stall 
forces for effective, sustained maneuvering in the high angle of attack, “post-stall” flight 
regime, or supermaneuverability, for next-generation fighter and attack aircraft [Acf. 3]. 
Historically, attempts to analvze such coinplex, unsteady behavior relied heavily on 
empirical data, obtained from often laborious and time-consunung tests. With the ad- 
vent Of the supercomputer, however, computation of actual viscous flow fields about 
nioderately complex computational models can now be numerically achieved in a matter 
of minutes through solution of the Revnolds-averaged Navier-Stokes equations [Ref 
4]. These solutions, in and of themselves, however, are not sufficient to proniote insight 
to the mechanics and physics involved in such flows. This additional requirement, for 
effective visual portrayal of the flow field, is satisfied by application of high periormance 
interactive computer graphics workstauons and associated software. 

The long range goal of the Computational Fluid Dynamics field is the development 
of a thoroughly verified computer code for unsteady aerodynamics which, among a 
invriad of other applications, will provide future aircraft designers with the opportunity 
to derive full advantage from application of the dvnamic stall phenomenon. [9o thus end. 
the Fluid Mechanics Laboratorv (F ML) of the NASA-Aines Research Center (ARC), 


in conjunction with the Naval Postgraduate School Department of Aeronautics aid 


Astronautics has planned an assortment of parallel and complementary oscillating airfoil 
windtunnel experiments and computer simulations. 

The current studv utilizes existing dynamic graphics packages to generate real-time 
projections of flow fields about a NACA-OO12 airfoil oscillating in pitch and experienc- 
ing deep dynanuc stall, as provided by a Navier-Stokes solver developed bv L. N. Sankar 
(Refs. 5.6]. A modified version of the Sankar code was submutted via the FML front end 
VAX, to the NASA-ARC Crav X-MP 48 computer, which output flow field solutions 
at specified intervals throughout the oscillatory cycle. From this data, graphics files were 
generated which, in turn, were submitted to the Graphics Animation Systern (GAS) 
software as developed by Sterling Software under contract to NASA-ARC, and run on 
an IRIS 3000-series graphics workstation. (Final output, for demonstrative purposes, 
was then transferred to video tape.) Thus, the results of the study were two-fold. First, 
procedures by which interactive computer graphics could be efficiently (in terins of both 
coniputer-time and man-hours) and effectively (in terms of information display) incor- 
porated as an analysis and verification tool for future fF ML studies, were developed, 
tested and refined. Secondly, the resultant graphics were utilized as a tool for the on- 


going verification of the Sankar code. 


г.) 


Il. DESCRIPTION OF THE SANKAR CODE 


Simulation of complex phenomena occurring in real fluid flows requires accurate 
ноль to the full Navier-Stokes equations. The Sankar Navier-Stokes solver, devel- 
oped for Blade-Vortex Interaction (BVT) studies, solves the two-dimensional, unsteady, 
compressible Euler and Navier-Stokes equations in strong conservation form, utilizing 
an alternating direction, implicit method as the time marching algorithm. A body-fitted 
cud with clustering in the normal dnection is utilized to discretize the flow field. 
lurbulent shear stresses are simulated with a two-layer algebraic eddv-viscosity model, 
with modifications as described later. The code may thus be utilized for solution of 


steady or unsteady, myviscid or viscous and laminar or turbulent flows. 


A. GOVERNING EQUATIONS 
[n Cartesian coordinates, the 2-D, unsteady, compressible Navier-Stokes equations 


in strong conservative form may be written as 


— 
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5g +6,E +6,F = Re '(5,R + 6,S) (2.1) 
where g, B Е. R and Š are four element vectors with entries corresponding, in order, to 
the equations of continuity, momentum (x- and y-) and energy. If non-dimensionalized 
such that all values of length, velocity (u and v), density ( о ) and total energy per unit 
volume (e) are normalized with respect to section chord length (c), free stream speed of 
Сора у, Пее stream density ( 9.) and p_a_, respectively, these vectors may be Writ- 
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The Reynolds number (Re), pressure (p), speed of sound (a) and stress terms ( Tm) are 


defined as follows. 
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Assuming Stokes’ hypothesis to be valıd, the constant 4 ıs defined as -2/3 a. Since only 
airflow 1s considered, the ratio of specific heats ( y ) 1s defined as 1.4 and the Prandtl 
number (Pr), as one. 

Neglecting viscous terms results in the Euler equations, as the right-hand side of 
equation (2.1) goes to zero. Use of the strong conservative form (Le, the солонаи; 
entrv in E, д(ри)/дх, is solved in its present form, vice the more mathematically correct 
form of pdu/dx + udp/ox) allows the identical conservation of physical quantities (mass, 


momentum and energy) when finite difference schemes are applied. 


B. TRANSFORMED GOVERNING EQUATIONS 

The flow field region must be discretized into a transformed, finite difference mesh, 
or computational plane, in order to allow numerical solution of the governing equations. 
The most efficient grids are rectangular in shape with regular spacing or spacing gradi- 
ents. They must, however, correspond to a flow field grid which provides high resolution 
(minimal spacing) in regions where gradients are large, particularly in the boundary layer 
and about the leading edge. In response to the coordinate transformation involved in 
the grid generation process, the governing equations, too, must be transformed. Dv de- 
fining £ and y as functions of the cartesian coordinates (time is not transformed), this is 
simply a 2-D mapping procedure accomplished by application of the transformation 


Jacobian, 
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the governing equation becomes 
ôq + ôE +ô, F= Re (5,8 + 5,5) (210) 
where 
q=qlJ (2.10) 
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C. GRID GENERATION 

Grid generation for discretization of the flow field region ах be accomplished « 
conformal mapping. algebraic methods or numerically solving a set of partial differential 
equations. The original Sankar code utilizes either of the latter two methods, while the 
version currently veetorized for use in this studv uulizes onlv the the algebraic method. 
which results in a body-fitted. sheared parabolic coordinate svstem or C-grid. Utilization 
of such bodv-fitted grids allows synehionous rotation ol the entire grid and aifoil sec- 
поп for dinanuc cases, Using simple tiigonometric relations. This coordinate system 
satısfies the generul giid requirements for sinootliness throughout and fine spacing in 
regions where high gradients exist, such as the boundary layer or leading edge. 

Normalized geometric airfoil shape data, in Cartesian coordinates. is input in table 
format to the code, which utilizes an interpolative procedure to compute additional 
points and smooth the surface. The trailing edge region is modelled as a vortex sheet 
shape, or “cut”, which smoothly leaves the airfoil, tangent to the mean cainber line at 
the trailing edge. Algebraic manipulation of the section surface and cut allows grid 
generation in the transformed ¢ - y plane. Figure 1 shows the resultant grid when 
mapped back to Cartesian coordinates. 

Uniform spacing in the transformed plane’s £ -direction results in fine spacing at the 
leading edge in the real plane, but coarse spacing in the trailing edge region. due to 
uniform increments across the axis (as opposed to the standard C-grid which results in 
a region of finer spacing at the trailing edge). 5 -spacing in both planes increases. ap- 
proximately exponentially, in directions normal to the airfoil surface, with the magnitude 
of initial spacing being user defined. Though easv to generate, this grid’s effectiveness 
in the analysis of certain flow cases has proven somewhat limited, due to its ccarseness 


in the trailing edge region [Ref. 6, p.44}]. 


D. APPLICATION OF BOUNDARY CONDITIONS 

Consideration of an airfoil impulsively started from rest in a fluid with uniform 
properties throughout, provides the initial conditions for solution of the parabolic Euler 
and Navier-Stokes equations (Eq. 2.9) in the computational piane. [n addition. bound- 
ary conditions and artificial dissipation terms are required in order to achieve accurate 
solutions. 

Three boundaries exist Which are of interest. the surface, the far-field boundary (grid 
linut) and the trailing edge cut. Boundary conditions at the surface are driven by the 


no-slip condition which, in response to viscous effects, requires the fluid velocity at the 
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Figure 1. Algebraic C-Grid in Cartesian Coordinates 


boundary to equal the boundary velocity. Thus, in this reference frame, u and v are set 
to zero on the solid surface. Since the surface ts considered adiabatic, both de/dy and 
Oe/O > are set to zero and, likewise, pressure distributions are determined from ópioj and 
cp/oz equal to zero at the solid boundary. (This is equivalent to neglecting stress con- 
tributions in the momentum equations, Wluch 15 acceptable when dealing with high 
Revnolds numbers flows.) Boundary conditions at the cut are driven by the necessity 
to avoid discontinuities in the “continuous” portion of the flow field. Since the gud 15 
extremely dense in the y -direction in this region, averages of the values of ио 
nearest mterior pots are assigned to pomts along the cut, allowing a smooth transition. 
Since the grid cannot econonicallv be generated large enough to avoid disturbance уе- 
locities at the far-field boundary, boundary conditions must account for the presence of 
the airfoil. Linear small disturbance theory 1s apphed to deternune perturbation veloci- 
ties at points along the boundary. which are then added to free stream conditions. 
Downstream boundaries are treated such that entropy changes can be convected out of 
the computational domain [Ref. 6. p. 29], allowing shocks and boundary-laver generated 
vorticity to pass through the grid. 

It has been found that spatial derivatives are sensitive to the decoupling of odd and 
even points which necessarilv occurs in central difference schemes. Thus resuits in the 
generation of high frequenev errors in regions of large pressure gradients, such as are 
present in shocks or about stagnation points. Due to the high Revnolds numbers en- 
countered in these flows, dissipation provided by the viscous terms are not enough to 
eliminate the errors, necessitating the addition of two artificial dissipation terms, em- 
bedded in the numerical schemes. An explicit artificial viscosity term 1s input by the 


user, With a proportional implicit term assigned by the code. 


Е. TURBULENCE MODELLING 

Since during the derivation of the Navier-Stokes equations. no assumptions ure 
made regarding flow type, these equations are instantaneously valid for both lanunar 
and turbulent flows. The large range of time and spatial scales encountered in turbulent 
flows, however, makes solution of the instantaneous Navier-Stokes equations virtually 
impossible, at this time. As a result, the equations are time- or enseinble-averaged. Use 
of such equations, in response to turbulence, results in additional turbulent shear stress 
terms or the Reynolds stresses (named for Oswalde Reynolds, who initiated turbulence 
studies in the 1880's), which effectively are increases 1n shear stresses due to turbulent 


motion. The difficulties associated with the existence of additional unknowis without 


any additional equations is described as the closure problem and is dealt with through 
the use of turbulence models. These are based on empirical knowledge of turbulent 
flows and thus, provided solutions will be approximate and require some confirmation 
from experiment [Ref. 7]. 

Ihe most common turbulence modelling method involves nuxing lengths. as devel- 
oped bv Prandtl. Assunung that turbulent fluctuations are essentially the result of ve- 
locitv perturbations between adjacent streamlines, and that all fluctuating velocity 
components at a given point are of the same magnitude, it can be shown that turbulent 
Ве Viscosity { 4, ) 18 proportional to the magnitude of the local vorticity ( w ) [IRef. 
8, р. 387]. In Prandtl’s mixing length model, the proportional term is the square of a 
characteristic length, related to fluid turbulence intensitv. Prandtl suggested this char- 


m@ereristic length ( / ) be treated as 
Pec Y 


where y 1s the normal distance from the fluid boundary and k is empirically obtained. 
Thus, the kev to obtaining accurate solutions when utilizing models of this type lies in 
the mixing length expression. | 

The Baldwin-Lomax two-laver eddy viscosity model, based on Cebeci's two-laver 
model, is presently used in the Sankar code. This model divides the boundarv laver into 
two regions, an inner laver and an outer layer, with separate methods for determining 
eddv viscosity used in each. The boundarv for the two lavers is defined as that point 
where eddy viscosities produced by the two methods match. The inner laver utilizes a 
mixing length model where eddv viscosity starts from zero at the wall and is defined bv 


the following. 
(Up inner = Ра! (2.12) 
The mixing length term 15 defined by 
= кур (215) 


where y is the normal distance from the wall, x is the Von Karman constant and D is the 


Van Driest damping function, given by 
р = [1 –ехр(–у ЈА“)] ເຊື 


where 4* is an empirical constant. The outer layer eddy viscosity model is defined by 


ແ... a KC OP ad) (2213) 
where K and C, are constants. The Kk  :notf intermittencvy function, given bv 
Um od 
F ie y) ss [! ປ С) ] (2.16) 


ensures that eddy viscosity approaches zero as the edge of the boundary laver is ap- 


proached and the flow assumes external characteristics. C,,, 18 a constant. F,,, 18 a 
function dciined the following relation. 
i - 2 
Fake annu а а n) a (2.17) 


Uais the magnitude of the velocity profile's velocity range. Fu. Is the maximum value 


provided by the following. 
Fy) =y|@|D (2.18) 
Ymax IS the value of y corresponding to F,,,4. Constants in the Sankar code are defined 
as follows: 
А = 26.0 
к = 0.4 
K = 0.0168 


ດ, > 1.6 
Ckteb = 0.3 
Cr = 0.25 


In order to account for turbulence outside the boundary layer, such as that which 
occurs in the dynamic stall process, a modified turbulence model is available which ef- 
fectively increases the outer model’s mixing length. In this case, F,,,, and },,,, are deter- 


nuned by the following. 


Ку 2 y^ |o|D (2.19) 
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The modified turbulence model was developed in response to early predictions of flow 
separation at high angles of attack and has been found to cause premature reattachment 
during the downstroke of dvnamic cycles. [ts use, therefore, is only recominended until 


frail onset (Rel. 8.p. 33). 


F. CODE STRUCTURE 

The code 1s structured such that the following user defined parameters are input 
from logical unit 05 (appended to the end of the code for Cray processing): grid size, step 
size (dt), artificial viscosity magnitude, mean angle of attack (o), oscillation magnitude 
(x,). angle for suspension of the modified turbulence model, reduced frequency (k). free 
stream Mach number (11), Reynolds number (Re), distance of the first 4 - contour from 
the airfoil wall, starting time, pitch and restart flags and airfoil geometry data. Added 
КОБИ tor the present study are tie number of time steps for the code to march on 
the present run, the number of steps between each plot (output interval) and the total 
number of steps and plots completed on all previous runs. Variables are then initialized, 
previous stored solution datasets are retrieved and read (to allow incorporation of all 
datasets, including those from the present run, into a single, combined file) and flow field 
Starting conditions are input. А series of subroutine calls then generate the grid 
(AIRFOL, SING, TABINT, WRAP), cluster grid pornts for viscous flows (CLUSTR. 
STRTCH), rotate the airfoil (ROTGRID) to its actual angle of attack and compute the 
initial metrics (METRIC). 

At this point, the code ıs fully initialized for commencement of the flow field sol- 
ution process, Which is conducted via an iterative loop. At each iteration, tume is first 
marched forward one time step, followed by computation of the time dependent values 
of angular velocity, angle of attack and step change of angle of attack, according to the 


relations: 
o — 2kM,, sinQk M.) (2.20) 


01-00-02) cosL 2k M (7 - dt) ] 


i=! 
da = 0, — 641 
The grid is then rotated, in the physical plane, by applying the the following relations 


at each erid point. 
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x = x cos(do) — y sin(a) (Zei) 
y=ycos(da) — x sin(dsavepha) 


Following recomputation of the metrics, the solution 1s computed by subroutine SLPS, 
the ADI-algorithm, which calls DISSIP, for computation of the explicit dissipation 
terms, STRESS and RESI, for computation of the inviscid and viscous terms, respec- 
tively, AMATI and MATRIXI, for computation of the Jacobian and tts inverse in the 
¢ -direction and AMAT2 and MATRIX2, for computation of the Jacobian and its ш- 
verse in the 7 -direction. Subroutine STRESS also calls subroutine EDDY tne ит 
lence model, for computation of the viscosity coefficient. WALLBC enforces the wall 
boundary conditions and finally, solution files, as discribed later, are generated. 

Prior to resumption of the loop, performance coefficients are generated by subrou- 
tines LOAD and CPPLOT. Surface pressure coefficients are obtained from the relation: 


ວ Po — Ps (2.22) 


a) 


where p, is the surface pressure and ບູ 15 the free stream pressure. Skin friction coefh- 


cients are a function of the wall shear stress (7,) according to the relations: 


E 
— À (2.23) 
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C = 
ty = ве МИ 
W = u, — vx ë J(u,x, + v,y;) 


The aerodynamic loads of lift, drag and moment about the quarter-chord, are computed 


bv the following relations. 


C; 2 C, cos(a) — C, sin(a) 
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С, = С, sin(x) — C, cos(a) 


Cn = [GLY зада + (x = хала] 
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C, and C, are the normal and tangential forces obtained from summing surface pressure 


and skin friction forces according to the following. 
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G. CODE MODIFICATIONS : 
ie 1ollowing modifications were made to the code, resident on the FML VAX. 
during the course of this study. 


1. Solution output conimands are placed within the loop in order to provide interval 
output. 


2. Read commands are placed at the beginning of the main program in order to allow 
storage of all solutions (previous restarts and the current run) in a single combined 
file. 


3. The airfoil section is rotated to the initial angle of attack, vice rotation of the free- 
stream direction. 


4. Restarts may be made from Plot3D (output) format as well as vector format. 


Additional Plot3D files (surface pressure and skin friction line plots) are output 
omesuorouume CPPEOT. 


6. Grid dimensions are made true variables. 


Eee seminputs are reformatted for ease of use. 


HI. DYNAMIC GRAPHICS GENERATION 


A. INTERACTIVE COMPUTER GRAPHICS 

Due to the complexities and tune dependent phenomenon associated with unsteady 
flows, attainment of acceptable clarity in flow field solutions requires complete descrip- 
tive information across fine grids, at a large number of minute time intervals. Thus, 
analysis and visualization of data generated by codes such as Sankar’s, is difficult due 
to both the volume of information provided and its temporal nature. Requirements 
identified for effective flow field visualization include maxinuzation of the bandwidth of 
information transfer, such that it closely matches the capabilities of the human eve, 
maximization of the quality of graphical displays, by enhancing key features and sup- 
pressing others, and maximization of the controllability of information [Ref. 9]. Addi- 
tional advantages in information transfer can be achieved through the use of redundant 
coding and structured displavs [Ref. 10]. [n response to these requirements, interactive 
computer graphics workstations have been evolved to complement the super-computer. 
Workstation capabilities, in terms of geometrical transformation and screen update dis- 
patch have been utilized by programmers to produce effective representations of flow 
field motion, through synchronization of coordinated data sets. Solution clarity 15 en- 
hanced by the high degree of spatial resolution and large range of colors afforded by the 
workstation displays. The interactive capabilities which currently exist for semi-coimplex 
flows serve to improve visual cues for display of three-dimensional data sets and allow 


immediate access to regions of interest. 


B. HARDWARE 

The NASA-ARC Fluid Mechanics Laboratory is currentlv equipped with an 
IRIS-3000 series workstation configured with 4Mbytes of display list memory and 
equipped with z-clipping and z-buffering hardware and a multiple key mouse. The [RIS 
display processing unit's refresh memory is arranged as a two-dimensional array (768 x 
1024), with row-column positions matching display pixel x-y coordinates. 24 bits are 
reserved for each pixel at eight bits (256 intensities) per color, resulting in a range of 
16,777,216 possible colors. When displaving dvnamic graphics, the refresh memory 15 
divided, in order to meet user demands, since flvback time (about 1.3 4 sec) 1s too short 
for complete picture update. This double buffering mode utilizes half the memory for 


refreshing the screen displavs, while the other half is dynamically updated. This ensures 


smooth motion on the display, but divides the number of bitplanes per pixel, in half. 
Thus, the number of colors available drops to 4096. In order to avoid a similar decrease 
in the range of colors available, color maps are utilized. In this case, pixel values in the 
refresh miemory are routed to an index or color look-up table (with each entry composed 
of 24 predefined bits, which define the color), vice direct routing to the intensity digital- 
to-analog converter. The IRIS transformation rate, via a single, composed transforına- 
tion niatrix in homogeneous coordinates, 15 80,000 coordinates per second, with points 
and lines generated at 3 million pixels per second, or 40,000 inches per second. Polvgons 
can be filled (flat shading) at the rate of 44 mullion pixels per second, or 70,000 square 
mches per second. These capabilities allow generation of most displays at the rate of 
approximately 10 screens per second, which, under ordinary conditions, 1s greater than 
the fusion frequency of the human eve. Interactive transformations are ordinarily ac- 


complished via the mouse. 


C. PLOT3D SOFTWARE 

Solutions provided by the Sankar code consist of multiple binary files. specifically 
formatted for input mto the Plot3D graphics software, authored by Pieter Buning for 
NASA. Two filetypes exist, which are labeled as XYZ or grid tiles and Q or flow 
quantity files. Both tvpes are three-dimensional matrices, with placement of data within 
the matrices corresponding to the row-column addresses of mesh points. XYZ-files have 
a depth of two in the two-dimensional case and provide cartesian coordinate definitions 
of each mesh point. Q-files have a depth of four in the two-dimensional case and pro- 
vide computed flow quantities at each mesh point, corresponding to the q-vector of 
equation 2.1. Headers for each file list free-stream Mach, Reynolds number, angle of 
attack and time. 

Additional flow field quantities at each mesh point are computed according to the 
following relations. [Ref. 11] 
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Based on these values, the Plot3D software will produce plots of scalar functions 
(densitv, pressure, temperature, enthalpy, energy, velocity, entropy, ınomentun and 
shocks) suitable for contour plots, vector functions (velocity, vorticity, momentum and 
pressure gradient), particle trace functions (particle traces and vortex lines), shock wave 
locations and grid functions (computational meshes and walls). 

Particle traces are generated using trilinear interpolation of values of the vector 
function inside a computational cell and second-order Runge-Kutta steps to advance the 
particle in space. Five steps are required for each cell. The algorithin for computing 
shocks computes the Mach number component in the direction of the local pressure 
gradient. Locations where this value decreases through 1.0 1s plotted as a shock. Two- 
dimensional stream functions are calculated by integrating the mass flow across a coor- 


dinate line. 


Output from Plot3D 1s available in a variety of formats and is a function of user- 
defined attributes. For dynamic plots, output is in the form of graphics files. formatted 
for further manipulation by animation soltware. Device independent plotting (DIP) is 
enabled through use of the ARCGraph binarv file format which allows data to be ma- 
nipulated bv a collection of function libraries and utilities developed by the Advanced 
Computer Research Center at NASA-ARC. Parameter data, along with graphics prim- 
itives (device independent op-codes) are contained in these files, allowing data manipu- 
lation by the Crav, IRIS and attached VAXes. Plot3D generation of graphics files 
requires that several attributes be defined. In order to reduce user-tasking, Plot3D will 
initiallv search for a filename-specific initialization communications file, which contains 
all required inputs. Input files must be read singularly, which makes necessarv the sep- 
aration of those combined plotting files received from tle Cray. Output graphics files, 
however, are once again stored in combined files in order to speed further processing. 
User-defined attributes include axis scales and extreme values, function, number of 
contours, line colors and types and wall definitions. 

Subroutine CPPLOT, within the Sankar code, provides grid (XYZ) files which con- 
tain plotting information for surface pressure coefficient aud skin friction coefficient line 
plots, corresponding to each output flow field solution file. These line plots are scaled 
and translated after separation from the combined files and plotted utilizing the grid 
function. Likewise, an angle of attack pointer plot, utilizing a sine wave format (gener- 
ated external to the code), 1s also provided. Q-files containing dummy variables (not 
utilized for grid plot generation) are also provided to complement the line plot 
XYZ-files, as required by Plot3D. 


D. GRAPHICS ANIMATION SYSTEM SOFTWARE 

The final step in the visualization process 1s the animation or synchronization of 
coordinated graphics datasets. The Graphics Animation Svstem (GAS) is a software 
package developed by Sterling Software under contract to NASA and provided to edu- 
cational institutions throughout the country. [t is device specific, running onlv on IRIS 
workstations, and is written in the C programming language under the UNIX operating 
system. Graphics files, in ARCGraph format, are read and stored in the IRIS’ display 
list memory as objects by the menu-driven program and assembled according to se- 
quence file instructions (stored transformation matrices). In order to smooth motion 
associated with geometric transformations, up to 30 linearly interpolated matrices are 


provided bv the program, between user-identified kevframes. Titles, legends and object 
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attributes are also available. Ultimately. flow field representations may be transferred 
fom eo tape or aon) (iim, Interactive viewing (ordinanly availabie in real-time) is 
ПОШ ШЕ Ма а шоо ми шетел "view data" window. 

A complete description of GAS and associated video hardware may be found in 
Relerence la. 


IV. RESULTS AND DISCUSSION 


A. PROCEDURAL CONSIDERATIONS 

Integration of the FML IRIS with related CFD software and the Cray X-MDP/48 
was tested bv plotting the flow field about an airfoil experiencing deep dvnanuc stall, as 
provided by the Sankar code. The advantages of plotting dynamic stall are two-fold. 
First, it is an extremely complex and time-dependent phenomenon which 15 difficult to 
visualize from discrete data sets. Thus, the accuracy of solutions for thus type flow field 
can only be determined by consideration of the conplete flow field in both time and 
space. Dynamic graphics are therefore ideally suited for its study. Secondly, the FML 
is currently heavily involved in studies of dynamic stall, which include verification of the 
Sankar code. Cray output 1s thus used to full advantage. 

In order to provide smooth animations, an extremely large number of plotting sets, 
provided at equal time intervals, must be generated by the code during the oscillatory 
сусје. Since access and transfer of this data 1s required several times during the graphics 
generation process, it was found that storage of datasets in combined files provided the 
greatest efficiency, especially in transfers between the Cray and IRIS. Software resident 
on the [RIS is then employed to separate the combined file into individual datasets. 
(Embedded in these programs are schemes to scale data, as necessary to control their 
ultimate on-screen positions, allowing, for example, placement of line plots in a specific 
corner of the display or plotting, for comparison, two flow fields side-by-side.) Individ- 
ual datasets are automatically assigned names which correspe::d to those contained іп 
Plot3D initiation files, also resident on the IRIS. These files currently contain com- 
mands for the processing of 50 datasets. Thus, the combined file separation and 
ARCgraph file production (Plot3X) process must be completed at intervals. The use of 
repetitive file names during this process, however, increases automation to such an ex- 
tent that only a limited number of user inputs are required to complete the process. 
ARCeraph files as provided by Plot3D are again in combined file format, allowing easy 


(a single user command) input into animation software. 


B. THE DYNAMIC STALL PROCESS 
A prerequisite to review of any code-provided flow field solution is consideration of 
available empirical information. The following dynamic stall characteristics have been 


observed for both harmonically oscillating airfoils (Ref. 13] and airfoils undergoing 
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monotonic angle of attack increases (Refs. 14, 15]. First, the airfoil stalls at an angle of 
attack which exceeds the static stall angle. Second, the corresponding normal forces and 
pitching moments exceed those attainable in the static case and, finally, a rapid change 
in pitching moment magnitude, termed as “moment stall” bv Harris and Pruyn [Ref. 16] 
occurs several degrees in azimuth prior to the normal force decrease, or “lift stall”, vice 
simultaneous occurence as in the static case. The process which results in these char- 
acteristics can be broken down into a series of chronological events, as depicted in Fig- 
Meme irom Reference (3. While the specific characteristics for a given airfoil are a 
function of free-stream Mach number, Revnolds number, reduced frequency, mean angle 
of attack and oscillation amplitude, empirical data obtained from a NACA-0012 airfoil 
atk = .15, Re = 2.5x 10° and a = 15’ — 10° sin(w2) 1s provided in the following discussion 
to illustrate the dynamic stall process. 

At point (a) of Figure 2, the static stall angle is exceeded without any detectable 
change in flow over the airfoil. The boundary layer remains thin with no evidence of 
flow reversal at the surface. At point (b), («a 2 19 – 20°) flow reversal appears at the 
surface. Over the majority of the airfoil, the boundary laver remains thin and attached. 
ИСЕ ear POON Ol the airfoil, however, the boundary layer thickens as the flow 
gradually decreases to zero velocity. At point (c), large eddies appear in the boundary 
laver. By point (d), flow reversal has spread over much of the chord, from the trailing 
edge to x'c = 0.3, With as much as 50% of the airfoil experiencing flow reversal, no 
changes in С, or C, are detected. At point (e), ( «x 2 23.4 ) a vortex forms. The 
boundary layer at the front of the airfoil abruptly breaks down (simultaneously from x.c 
= 0.0 to 0.3) and the vortex moves downstream at approximately 35 - 45% of the free- 
stream Velocity. At point (f) the lift-curve slope increases beyond the 2zo limit of 
quasi-static flows. By point (g), the pressure distribution is altered sufficiently to 
produce a noticeable divergence in the pitching moment (moment stall) but 1s accom- 
panied by a continued increase in lift. Maximum lift occurs at approximately mid-chord, 
followed by a sharp decrease (point (h)), which identifies lift stall. (Boundary laver sep- 
aration has occurred to such an extent that continued increases in angle of attack will 
not result in continued increases in lift.) At point (i) ( « — 24.95 ? ), maximum negative 
moment occurs. The vortex moves off the trailing edge at « = 24.8 * (downstroke) and 
C, and C, move toward static levels. At point (j), the airfoil is fully stalled. At point (k), 
boundary laver flow begins to reattach, progressing from the leading edge at approxi- 


mately 25 - 35% free-stream velocity and is complete by a = 7 ° (downstroke). At point 
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INCIDENCE, a, deg (1) RETURN TO UNSTALLEO VALUES 


Лоше 2. The Dynamic Stall Process (from Carr, Ref. 13) 


(1), after some hvsteresis ( « = 6° , upstroke) the separated region has fully closed and 


unstalled values of C, and C, are reestablished. 


C. PROCEDURAL VERIFICATION 
The case chosen for procedural verification closely corresponds to the Carr data of 
the previous section and exactlv matches conditions previouslv run when the code was 
imually vectorized for use on the АКС Cray (Ref. 17]. There are several advantages 
associated with this duplication. First. it allowed storage of a complete record of this 
case (data saved at 300 intervals, vice 12) for future access. Second, it provided a simple 
means bv which to determine if code modifications were correctly incorporated. Finally, 
It provided a baseline, in conjunction with empirical data, against which future sensitiv- 
itv studies could be compared. 
Inputs for this case were as follows: 
не = 35x 105 
И 293 
Ot = .005 sec. 
я -тш = .00005 


с) = .131 


Di 


|| 
Сл 


Artificial Viscosity 


Grid size: 157 x 40 


Oscillatory Cvcle: a = 15° — 10° cos(we) 


The original Baldwin-Lomax turbulence model was used throughout the oscillatory 
cvcle. Figures 3 through 7 are coefhicient of pressure contour plots of the predicted flow 
field. (Dynamic plots are similar to these Plot3D-provided plots, but do not contain 
explicit quantitative information.) As previously noted [Ref. 17, p. 48], under these 
conditions, moment coefficients are consistently underpredicted as compared to exper- 
imental data. Drag and lift coefficients closely match experimental data below approx1- 
mately 15 and 18 degrees angle of attack, respectively. Early prediction of flow 
separation, however, forces inaccuracies in quantitative solutions bevond these values 
and underprediction of the maximum lift coefficient obtainable. These findings highhght 
the rationale behind the modified turbulence model. Information provided bv the dv- 
namic plots, shows general qualitative agreement with experimental data, during the 
upstroke. A vortex forms near the leading edge of the airfoil section and progresses 


down its upper surface, graduallv increasing in size. As the downstroke begins, however, 
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rather than shedding from the trailing edge, the vortex stagnates at approximately the 
65% chord point and gradually dissipates. Line plots of surface pressure coefficients 
more accurately portray (as compared to traditional carpet plots) the turbulent effects 
associated with the existence of the vortex. During the vortex dissipation phase, these 
line plots show a slight (gradually decreasing) pressure rise in the vicinity of the vortex. 
The subtle differences between these plots and those which would have been obtained 
had the vortex actually shed, are lost when this information is integrated for lift. drag 
and, to a lesser extent, moment coefficient data. Thus, quantitative data may not be 


good indicators of flow field solution accuracy. 


D. SENSITIVITY ANALYSIS 
The following code parameters which should effect dissipation of the vortex have 


been identified. 


|. Artificial viscosity magnitude. 


(~) 


Turbulence model parameters. 


Grid resolution. 


us 


A limited sensitivitv analvsis into the effects of grid resolution changes in the trailing 
edge region was conducted in order to provide direction for follow-on studies. During 
the initial phase of this analysis, only qualitative results will be considered, for those 
reasons mentioned above. Dynamic graphics, therefore, are ideally suited for this task, 
especially when data scaling capabilities are utilized in order to plot comparative data- 
sets, side by side. 

The grid ts relatively coarse in the region in which the vortex becomes stationary and 
dissipates, as compared to those regions in which it behaves as expected. In thie first 
sensitivity case, therefore, grid dimensions were held constant and the distance of the 
first constant- y line from the airfoil surface was increased to .001 , resulting in a finer 
grid in the trailing edge region. Figure § shows the effects of this change on the vortex. 
Within this grid, the vortex once again becomes stationary, but effectively splits, such 
that a secondary Vortex is shed while the original vortex dissipates. While this does not 
correspond to any process encountered in actual flows, it does indicate that final sol- 
utions are sensitive to this parameter. 

Changes in y -spacing have two effects on the flow field solution. Grid spacing in 
the trailing edge region becomes more fine at the expense of grid spacing in the boundary 


laver, which becomes more coarse. The extent of influence of the boundary layer sol- 
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ution on the vortex. or far field solution, is not yet clear. [t is known. however, that 
boundary laver solutions are extremely sensitive to y -Spacing parameter changes, re- 
quiring an initial value of .QOUOS for accurate solutions (Ref 18|. Further study, there- 
fore, is required in order to determine whether vortex behavior changes observed in the 
previous case were in response to the altered boundary layer solution (which would re- 
commend investigation of turbulence model parameters) or the altered grid. In order to 
facilitate this study, the existing grid would require enlargement (from 161 x 41). Un- 
fortunately, this results in such a large memory allocation from the Cray, that it is pro- 
hibitivelv inefficient. Replacement of the entire grid generation process is therefore 


currently under consideration for this code. 


E. CONCLUSIONS 
The current study has resulted in completion of the following items. 


[. Full procedural documentation (Appendix A) has been developed for efficient, 
code-independent, two-dimensional, dvnanuc graphics production on FML and 
NPS IRIS workstations. Data must be in Plot3D format and may be provided bv 
either code or experiment. A method for simultaneous viewing of multiple datasets 
has been incorporated. 


2. The Sankar Navier-Stokes solver (Appendix B) has been modified co allow dynamic 
graphical presentation of its flow field solutions. 


3. Dvnamic graphics applications in the studv of complex flows were illustrated bv 
means of an introductory sensitivity analysis, using the Sankar code. This limited 
analysis identified either grid resolution in the trailing edge region or the boundary 
laver solution as parameters to Which vortex behavior is strongly dependent. (Thus 
suggests that incorporation of a more sophisticated grid generation scheme would 
be worthwhile.) 

In order to derive full benefit from application of this technologv, existing and future 
empirical data should be stored in plotting forinat, allowing the availability of visual re- 
cords of multiple flow functions. Also, modification of additional Navier-Stokes solvers 
for dvnamic output will allow effective code-to-code and code-to-experiment compar- 


isons, all in an effort to eventually develop an accurate flow field predictor. 


3] 


APPENDIX A. PROCEDURES FOR GENERATING DYNAMIC 
GRAPHICS AT FML 


A. THE ANIMATION PROCESS 
Creation of dynamic graphics for flow field visualization at the NASA-Ames Re- 
search Center Fluid Dynamics Laboratory involves five basic processes: 


1. Navier-Stokes solver is submitted to the Crav A-MP,48 and dynamic output saved 
in the Central Storage Facility (CSF) in combined files. 


2. Combined files are converted to IRIS binary and transferred to an FMLIRISI ac- 
count. 


2 


Combined files are separated into individual plotting files on the IRIS and con- 
verted to combined graphics (.gra) files bv Plot3D. 


4. Graphics files are fed to the Graphics Animation System software (GAS) and dy- 
namic graphics are plotted. 


Сл 


Resultant files and plots are transferred to storage or video tape. 


Prior to beginning the animation process, the code must be modified to provide dy- 
namic output in Plot3D format. This involves: 


1. Placement of solution output commands within the iterative loop along with some 
interval flag. 


гә 


[f restarts are to be used, placement of read statements at the beginning of the 
program to allow storage of all provided solutions in a single combined file. 


3. Storage of any additional required plotting information (line plots) in Plot3D XYZ 
files for eventual display utilizing Plot3D’s grid function. 


Full documentation of Plot3D formats and functions may be found in Sterling Software 
technical note 15, A Guide to Generating Movies using Plot3D and GAS, by Greg Howe. 
The procedures outlined below refer specifically to a version of the Sankar Navier-Stokes 
solver, modified for dynamic output. They can easily, however, be generalized to apply 
to any code (or experiment) from which dynamic output 1s desired. While many meth- 
ods of completing these steps may exist, those outlined below have proven to be the 
most time efficient. 
i. Vax Submitted Cray Jobs 

A full dynamic cycle, utilizing a 161 x 41 grid, requires over 5000 seconds of 

Cray time, making job chaining necessary to obtain a reasonable priority. With 900 


second jobs, 24 to 48 hours are usually sufficient for the full computation to be com- 


pleted. Since significant differences exist between the first JCL from which Plot3D files 
are saved and subsesquent restart JCL's, two JCL's are maintained on identical copies 
ol Samkar's code. 

ОЛОИ lee ere Necesses 5 Of Chdv tape (vid Si GIN commands) for initial 
ОО ОМА бс е ОВО ог тпагах format. Combined Plot3D files are 
initialized on CSF and assigned PDN’‘s and ID's. Tape 05 input data (appended to end 
of code) is assigned in accordance with definitions included in MAIN portion of the 
ОО restart 15 Пош Ч патос data, ADA will not be a whole number. TSTART 
Mest be adjusted accordingly, ш order to have accurate time ADA matches.) 

RESTART.JCL. Таре US inputs must idenucally match those from 
met JC, with the exception of ISTART,CSTP,CPLT and possibly FORMAT. 
Peel. JCE-defined PDN ID solutions are accessed (read, stored and appended with 
additional sol.tions), saved (witli identical PDN ID's} and scrubbed (oldest versions 
deleted). The initial restart solution is accessed with library routine GETP3D. written 
puc owe which will read a combined file, skip a specified number of datasets 
me ool} and access the neXt dataset (or datasets, 1f NUSMDS is not cqual to one, the 
default). Thus, 1f 

Dookie + GUN | 


the correct solution will be provided for restarts. Subsequent restarts illustrate the ad- 
Mintase rol combined file name (PDN ITD) repetition, as they simply require the following 
feo lA JCL and Tape 05 updates: 

Peo sok IP 

2. Job chain commands. 
CS ly 
4. CPLT 


(,2 


Once the optimum number of plotting sets for a cycle 1s deternuned, plotting 
files for an AOA pointer on IRIS displays can be generated using Plot3D. (500 sets 
maximum, or modify dimension statements.) 

Line Plots. Any pointer or line plot information may be generated, either within 
the code or externally, and stored in XYZ files for further processing. Utilization of 
Plot3D function 0 and proper descriptors of walls will produce the desired plots. 
SINE.TXT is a program which generates plots for a dynamic AOA pointer in sine wave 


format. In this program, Tape 05 includes values for pointer scaling (axis length with 
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no scaling is 27 ) and positioning on the IRIS screen and initial pointer locatiou. For 
example. if the cycle starts from the mummers 


PNLOCI = ВОВ 


Hardcopy printouts mav also be obtained fioi these grid files by use of pro- 
erams such as CPPLOT.TXT. This program provides printouts of upper and lower 
surface pressure coefficient and skin friction values at x € locations, including Cp plot 
generation. 

2. Path to the IRIS 

SEND.JCL. Ctilizes GETP3D and SENDWAS to convert to [RUS pir ee 
transfer complete or partial combined files, or individual plotting sets, fron tlie CA Y 
to the IRIS. (Since I O's between the Cray and IRIS can get clogged ГЕ 
transfers are requested, the individual method is net recommended.) This max need to 
be completed in stages [0 а о ексвеф те ле ЕТТЕ 

3. Creating Graphics Files 

Combined files on the [RIS are separated into individual plotting sets by the 
programs GETSIN.F (also generates the dummy Q file, “ysin.iris”), GETCP.F (also 
scales and positions Cp or skin [friction on the [RIS screen) and СЕ ТХ Е (Tor DEN 
Q flow field datasets). These programs require some initial user irıputs prior to running. 
Output dataset names will automatically match Plot3D imtialization file names. Since 
the Plot3x run for each dataset type (AX, P and SIN) requires specilic arguments, 1nitlal- 
ization files for each data type exist separately in files AZINECOMFETN Gs ae 
SINI.COM. Since Plot3X searches for the fle name PLOTSDENELCOM F АНЕ 
must be renamed appropriately, utilizing the “mv” command. Enteriig the command 
“plot3x” wil then generate multiple graphics files which are stored in a сотошешие 
matching the initial Q file name encountered by Plot3D, (1.e., (001 ота) ІШПЕ 
then stored on an appropriate subdirectory or fed to GAS. 

4. Notes on GAS 

For movies, the maximum nuniber of objects per sequence is Írfiv. The 
“seq(n).seq” files are 50 object far-field solution files (seq(n).seq = objects п(1-30)), 
which use object 400 for titles. They can be fed to GAS (after object input) via the ALA 
[О window. For interactive viewing, the VIEW DATA window is the only usable op- 
tion. (No sequence files are required in view data.) Full GAS documentation is avail- 


able in the GAS User's Guide, available from Steiliug Software. 
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5. Storage 
Since [RIS memory space is severely limited at FML., resultant graphics files 
should be stored elsewhere, if not in immediate use. Downloading to ] 4” cassette tape 


pxeasibaacconplished (also a reconunended backup). as 1s traasfer to CSF. 


B. SUPPORT CODES 
il. INITIAL.JCL 


JOB, JN@FLYNAVY, T=90@.  ERIC PAGENKOPF , X4269. 
ACCOUNT . AC= , US= , UPW= . 

* . 

CFT ONA OFFS: 

.. 

e RESTART COMMANDS: 
CSF ACCESS, CSFACCT=,CSFPSWD=,DN=0LDSLN , PON=NSE509. ID=VALDES. 
ASSIGN,DN=OLDSLN,A=FTO7. 

• САМЕ DATA FROM PRESENT RUN IN COMBINED FILE: 
ASSIGN,DN=XYZ, А-ҒТ280. 

ASSIGN, DN=Q, A=FT21 . 

ASSIGN,DN=CPX,  A=FTSO. 

ASSIGN.ON=CFX,  AmFT6Q. 

* 


LDR .MAP=PART , SET=ZERO . 











CSF ,SAVE ,CSFACCT= , CSFPSWD= .DN テ XYZ , PDN=SNKROSX, ]D=FLYNAVY. 
CSF,SAVE,CSFACCT=e,CSFPSWD= , DN=Q, PDN=SNKROSQ, ID=FLYNAVY. 
CSF,SAVE,CSFACCT=,CSFPSWD= , DN=CP X , РОМ=5МККӨ5Р, ID=FLYNAVY. 
CSF , SAVE , CSFACCT=,CSF PSWD= ,DN=CFX, PDN=SNKROSF , ID=w=FLYNAVY. 


* . 


. ACCESS , ON=SENDVAX , PON=SENDVAX, ]D=STTROM , OWN=RFTRDM . 





s. SENDVAX,DNzXYZ , VONS'RAL"JIANPS password": : [JIANPS.PAGAN.DATA]XT1.DAT'. 

s. SENDVAX,DN-Q, VON 'RAL"JIANPS password" : : [JIANPS.PAGAN.DATA]QT1.0OAT'. 

s. ACCESS, DNSSENOWKS , PONZSENDWKS , 1D=STTRDM, OWN=RFTRDM. 

«. SENDWKS, DN-Q3X49, VON 'FMLIRIS1"jiaa possword"::"/u/jiaa/pagan/q. iris" ' , NOREC. 
e. SENDWKS , DNaCFX, VONz'FMLIRIS1"jiaa password"::"/u/jiaa/pagan/cf01.dat"' , NOREC. 
ມ ລພ <<- < < < ະນ — a o Fa mmawusasz===— og EEE i ABE ລສທ ເສສ > 


.. „ОВ СНАТМ1 НС COMMANDS: 
FETCH,DN=JOB2,TEXT='[PAGANJ]NSMULT.TXT; 2°. 
REWIND,DN=JOB2. 

SUBMIT,DN=JOB2. 


“. 
/EOF 
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ВЕКЕ. СТ. 


JOB, JN=FLYNAVY,T=988. ERIC PAGENKOPF ,X4269. 
ACCOUNT , AC= , US=, UPW= . 

a 

CFT,ON=A ,OFF=S. 

е. 

e. RESTART COMMANDS: 
e.CSF,ACCESS,CSFACCT=,CSFPSWD= , DN=OLDSLN, PDN=NSES89, ID=VALDES. 
CSF,ACCESS,CSFACCT=,CSFPSWD= ,DN=GETP3D, PDN=GETP3SD, ID=RFRICC. 
CSF,ACCESS,CSFACCT=,CSFPSWD=,DN=XYZ1, PDN=SNKRO5SX, [DF LYNAVY . 
CSF, ACCESS, ,CSFACCT=,CSFPSWD=,DN=Q1, PDN=SNKR@5Q, ID=F LYNAVT . 
CSF ,ACCESS.CSFACCT=,CSFPSWD= , DN=CPX1, PDN-SNKRO5SP, IDEFLYNAVY. 
CSF , ACCESS , CSF ACCT= ,CSFPSWD= , DN=CFX 1 。 PON=SNKROSF , ID=FLYNAVY, 
REWIND ,ON=Q1. 

GETP3D, I=Q1,0=0LDSLN,DATATYPE=Q,DSSKIP=49. 

REWIND,DN=OLDSLN. 

REWIND,DN=XYZ1. 

REWIND, ON=Q1. 

REWIND ,ON=CPX1. 

REWIND ,ON=CFX1. 

ASSIGN,DNzXYZ1, АжҒТ90. 

ASSIGN, ОМ№=01, A=FT91. 

ASSIGN,DNzCPX1,  Ax=FT92. 

ASSIGN,DN=CFX1,  AsFT93. 

ASSIGN, DON=OLDSLN,AmFTO7. 

е . 
e SAVE DATA FROM PRESENT RUN IN COMBINED FILE: 
ASSIGN, ON=XYZ, AmFT20. 

ASSIGN, ON=Q, AmFT21. 

ASSIGN, ON=CPX, АеҒТ50. 

ASSIGN,DN=CFX, AmFT60. 








.. 

LDR , MAP=PART , SET=ZERO . 

Ф, 

CSF ,SAVE ,CSFACCT=, CSFPSWD= , DN=XYZ , PDN=SNKROSX, ID=FLYNAVY . 
CSF, SAVE, ,CSFACCT=,CSFPSWD=,DN=Q, PDN=SNKROSQ, ID=FLYNAVY. 
CSF,SAVE,CSFACCT=,CSFPSWD= ,DN=CPX, PDN=SNKROSP, ID=FLYNAVY . 
CSF, SAVE, CSFACCT=,CSFPSWD= ,DN=CFX, PON=SNKROSF , ID=FLYNAVY. 
CSF , SCRUB,CSFACCT=,CSFPSWD=,PDN=SNKROSF, ID=FLYNAVY. 

CSF , SCRUB , CSFACCT= ,CSFPSWD= , PDN テ SNKR95P , ID=F LYNAVY . 

CSF ,SCRUB,CSFACCT=,CSFPSWD= ,PON=SNKRO5X, ID=FLYNAVY. 

CSF ,SCRUB,CSFACCT=,CSFPSWD= , PON=SNKROSQ, 1 D=F LYNAVY. 

e ACCESS , ON=SENDVAX , PON=SENDVAX , 1D=, OWN=. 

. SENDVAX , DNzXYZ, VDN='RAL"JIANPS passwoard":: [JIANPS.PAGAN.DATA]XT1.DAT'. 
.SENDVAX,DN-Q, VDNZ'RAL"JIANPS passward": :[JIANPS.PAGAN.DATA]QT1.DAT'. 

„ACCESS , ON=SENDWKS , PON=SENDWKS, 1 D=, OWN=. 

. SENDWKS , ON=03X49,VON=’FMLIRISI"j ica passward"::"/u/j iaa/pagan/q.iris”’ ,NOREC. 
.SENDWKS , DN=CFX,VON='FMLIRISI”"jiaa passward"::"/u/j iaa/pagan/cf@1.dat"’ ,NOREC. 








e. JOB CHAINING COMMANDS: 
FETCH, DN=JOB3, TEXT=' [PAGAN]JOB83.TXT’. 
REWIND, ON=JO83. 
SUBMI T , DN=JOB3 . 





в. 
/ЕОЕ 


у] 


3. SEND AGE 


JOB, JN=GETX,T=38. ERIC PAGENKOPF, X4269. 
ACCOUNT , AC= , US= , UPWz . 


CSF ACCESS .CSFACCT=,CSFPSWD=,DN=XYZ , PON=SNKR38X, I D=F LYNAVY . 
CSF ACCESS .CSFACCT= ,CSFPSWD= , DNeQ , PDN=SNKR38Q , ID=FLYNAVY . 


CSF . ACCESS .CSFACCT= .CSFPSwD= , DN=CETP3D , PDNeGETP3D . ID=RFRICC.。 
ACCESS. DN=SENDWKS, PDN=SENDWKS . ID= , OWN= . 


REWIND, DN=XYZ. 
REWIND .ON=Q. 


GETP3D. I=XYZ .O=XOUT 1 .DATATYPE=XYZ . DSSKIP=169 . NUMDS=156. 
GETP3D. l=Q, O=QOUT1, DATATYPE=Q, DSSK [P=10@, NUMDS=15@. 
SENDWKS, DN=XOUT1, VDN='FMLIRIS1"jiaa password": :"/u/j iaa/pagan/38x.iris"' ,NOREC. 
SENDWKS,DN-QOUT1,VDN-'FMLIRIS1"jiaa passward"::"/u/jiaa/pagan/38q. iris" ' , NOREC. 
« .FETCH,DN-JOB2, TEXT&' [PAGAN]SEND. JCL;2'. 
* REWIND, ON=JOB2. 

¢ . SUBMIT, DN=JOB2. 
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JOB , JN=F LYNAVY , T=30. 


4. SINE.IXT 


ACCOUNT , AC= , US= , UPW= . 


ERIC PAGENKOPF , X4269. 


ACCESS , DN=SENDWKS , PDN=SENDWKS , ID= , OWN=. 


CFT.ON=A.OFF=S. 


ASSIGN.DN=X, A=FT50. 
ASSIGN,DN=Q, A=FT60. 


LDR , MAP=PART , SET=ZERO. 


SENDWKS .DN=Q , VDN=' FMLIRIS1"jioo PASSWORD": :/u/j iaa/pagan/qsin.iris’,NOREC. 


SENDWKS , DN=X 


.VNDNe'FMLIRIS1"jioo PASSWORD"::/u/jioo/pogon/snkr05s' , NOREC. 


CSF,SAVE,CSFACCT=,CSFPSWD=,DN=X , PDN=SNKROSS, ID=FLYNAVY . 


/EOF 


PROGRAM POINTER 


One а е € € а а €9 а а 9 е е е а а а 


Caso 


C 
C 


Case 
C 


20 


THIS PROGRAM GENERATES PLOTSD FILES WHICH WILL PRODUCE AN 
AOA POINTER WHEN FUNCTION 0 
ON CSF FOR FUTURE ACCESS VIA PROGRAM GETSIN.JCL. A SINGLE Q 
FILE OF PROPER DIMENSION IS SENT DIRECTLY TO THE IRIS WITH 


IS SELECTED. XYZ FILES ARE SAVED 


FILENAME QSIN.IRIS. 


VARIABLES: 


NPLOTS: TOTAL NUMBER OF PLOTS IN ONE CYCLE. 


PNLOC 
PNLOC 
XSCAL 
YSCAL 


1: INITIAL INTEGER LOCATION OF POINTER 
: INTEGER LOCATION OF POINTER. 

E: AXIS SCALING FACTOR 

E: SINE WAVE SCALING FACTOR 


XPOSIT: X-POSITIONAL FACTOR (SCREEN LOCATION) 
YPOSIT: Y-POSITIONAL FACTOR (SCREEN LOCATION) 


INTEGER PN 


DIMENSION SINX(3,500),SINY(3,500),01(3,500),02(3,500) 


LOC.GDIM,PNLOCI 


DIMENSION Q3(3.500),04(3,500) 


DATA Q1,Q2 
INITIALIZE 


READ (5,1) 
READ (5,2) 
Р]=АТАМ(1. 
GDIM = 3 

PNLOC = PN 


DO 19 L=1, 


,Q3.Q4/6999* 1 ./ 


NPLOTS,PNLOCI , XSCALE, YSCALE,XPOSIT,YPOSIT 
jes. 


LOC] 


NPLOTS 


GENERATE PLOT3D POINTER FILE 


DO 20 N=1,NPLOTS+1 


XLOC = Мз2. еР] /МРЕОТ$ 


SINY(1,N) = (0.0)eYSCALE + YPOSIT 
SINX( 1.N) = (XLOC)*XSCALE + XPOSIT 
SINY(2,N) = (SIN(XLOC))*YSCALE + YPOSIT 
SINX(2,N) = (XLOC)=XSCALE + XPOSIT 

CONT INUE 


PLOC = PNLOC+2.+PI/NPLOTS 
SINY(3,1) = (@.@)*YSCALE + YPOSIT 
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SINX(3,1) = (PLOC)*XSCALE + XPOSIT 
SINY(3,2) = (SIN(PLOC))*YSCALE + YPOSIT 
SINX(3,2) = (PLOC)*XSCALE + XPOSIT 


(> 
WRITE(59) CDIM,NPLOTS 
WRITE(SO)((SINX(1,3), I=1,3),J=1,NPLOTS), 

+ ((SINY(I,J), I=1,3),J=1,NPLOTS) 

C 
PNLOC = PNLOC + 1 
IF(PNLOC . EQ .NPLOTS ) THEN 

PNLOC = PNLOC + 1 - NPLOTS 
ENDIF 

10 CONTINUE 

C 

Ceee GENERATE DUMMY Q FILE 

C 

МВ 1ТЕ( 60) СОМ, МРЕОТ$ 

WRITE(6@) GDIM,GDIM,GDIM,GDIM 

WRITE(6@) ((Q1(I.J), 1=1,3),J=1,NPLOTS), 
+ ((92(1,4/), I=1,3),J=1,NPLOTS), 
+ ((Q3(1,J), I=1,3),J=1,NPLOTS), 
+ ((04(1,J), I=1,3),J=1,NPLOTS) 

C 

1 FORMAT (1X) 

2 FORMAT (1Х,218,4Ғ8.4) 

END 
/EOF 
NPLOTS: PNLOCI: XSCALE: YSCALE: XPOSIT: YPOSIT: 
294 222 . 98 . ] -1.0 -.5 
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Di 


св ВОТ. ГАГ 


JOB , JN=SKYHAWK , T=30. ERIC PAGENKOPF ,X4269. 
ACCOUNT , AC= , US= , UPW= . 


е. 
CFT,ON=A ,OFF=S. 


CSF, ACCESS , CSFACCT=, CSFPSWD= , DN=GETP3D, PON@=GETP3D, ID=RFRICC. 
ACCESS , DN=SENDWKS , PDN=SENDWKS , IDzSTTRDM, OWNERFTRDM. 


CSF, ACCESS, CSFACCT=,CSFPSWD= ,DN=CPX1, PDN=SNKROSP, ID=FLYNAVY. 
CSF,ACCESS,CSFACCT=,CSFPSWD=,DN=CFX1, PDN=SNKROSF, ID=FLYNAVY. 
CSF,ACCESS,CSFACCT=,CSFPSWD= ,DN=0Q1, PDN=SNKROSO, ID=FLYNAVY . 


REWIND , DN=CPX1. 

REWIND ,DN=CFX1. 

REWIND ,ON=Q1. 

GETP3D, I=CPX1,0O=CPX,DATATYPE=XYZ,DSSK]P=49 ,NUMDS=2. 
GETP3D, I=CFX1,O=CFX,DATATYPE=XYZ ,DSSK1P=49 ,NUMDS=2. 
GETP3D,I=01, O=Q, DATATYPE=Q, DSSKIP=49 ,NUMDS=2. 
REWIND , DN=CPX 

REWIND .DN=CF xX. 

REWIND , DN=Q . 

А5510М,ОМеСРХ, А«ҒТ50. 

ASSIGN,DN=CFX,A=w=FT60. 

ASSIGN,DN=Q, Am=wFT70. 


9, 
LDR , MAP=PART , SET=ZERO . 
.. 
/EOF 
PROGRAM CPPLOT 


-еееегеееееееееееееееееееееееееегееееееееееееееегееееееееегеееееееееееееееее 


& 
е THIS PROGRAM READS STORED, TRUE VALUES OF XOC, CP AND CF, 
“ AND PRINTS THEM OUT TO TAPE 06. 
Ф 
ә VARIABLES: 
» NUMDS: NUMBER OF DATA SETS IN THE READ FILE 
e NPLOT: PLOT NO. OF FIRST DATA SET (DSSKIP + 1) 
Ф 
a TAPE ASSIGNMENTS: 
è TAPE 05: INPUT DATA 
a TAPE 06: OUTPUT DATA 
. ТАРЕ 50: TRUE CP VALUES 
è TAPE 60: TRUE CF VALUES 
è TAPE 70: Q FILES 
. 
の ooo の の の の の の の の の の の の の の の の の の かよ の の の ひひ のか の みみ の かよ の の の の ゆい の ゆ ゆ の の の の よ ゆか ゆ の の の ゆ の か お の かみ の の の の の の ゆ の の の の の の の の の ゆ ゆ の の の 
C 
DIMENSION CPX(3,49),CFX(3,49),CPY(3,49),CFY(3,49) 
DIMENSION Q1(161,41),02(161,41),03(161,41),04(161,41) 
DIMENSION KODE(4),LINE(90) 
DATA KODE/1H ,1H+,1HI,1He/ 
С 
Ceee READ INPUT DATA & INITIALIZE 
С 
КЕАО (5,1) 
READ (5,2) NUMDS,NPLOT 
DO 9 Iz1,90 
LINE(I) = KODE(1) 
9 CONT I NUE 
С 
Cees READ & SAVE TRUE VALUES 
С 


41 


* 
a 
a 
9 
a 
* 
® 
3 
ёа 
ж 
а 
з 
Ф 
$ 
o 
Ф 


1 


O Y Lin — O) — 


DO 10 L121,NUMDS 
ADS IMAX , KMAX 
READ (50 າ I=1,IMAX), 
+ ((CPY(1,J 
КЕАО (60) ІМАХ,КМАХ 


. 11, ТМАХ ), 


Ј=1,КМАХ), 
J=1 , KMAX ) 


КЕАО (60) S ‚ Imf, IMAX),J=1 е" 


+ ((СЕТ(1,Ј), 1=1, 1МАХ) ,4=1,КМАХ 
ສ. ЈМАХО , КМАХО 
READ( 79 AMINF ,ALFAD,REYREF , TIME 
КЕАО (70) :(Q1(I,J). Im1,IMAXQ), Jz1,KMAXQ), 
+ ((Q2(1,J), I=1,IMAXQ),J=1,KMAXQ), 
+ ((Q3(1,J), I=1.IMAXQ),J=1,KMAXQ), 
+ ((О4(1.Ј), I=1,IMAXQ) ,J=1 ,KMAXQ) 
CPO = (1. + .2 © AMINFee2) se 3.5 - 1. 
CP@ = СРД / ( .7 « AMINFee2) 
KO = 30. e CPO + 4.5 
DO 11 L=1,KMAX 
K = 30. e (CPO — CPY(3,L)) + 4.5 
КІ к 30. » (СРО -СРҮ(2,1)) % 4.5 
IF(K.LT.1) K = 1 
IF(K1.LT.1) K1 = 1 
IF(K.GT.90) K ж 90 
IF(K1 .GT. 90) КІ = 90 
LINE(K@) ж KODE(3) 
LINE(K) = KODE(2) 
LINE(K1) = KODE(4) 
IF(L.EQ.1)THEN 
WRITE(6,*«)' PLOT AOA TIME MACH REY NO. 


WRITE(6,3) NPLOT,ALFAD,TIME,AMINF,REYREF 


WRITE(6,*)' XOC CFL CFU CPL 


WRITE(6,4) CPX(1,L),CFY(2,L),CFY(3,L) 
ELSE 

WRITE(6,4) CPX(1,L),CFY(2.L),CFY(3,L) 
ENDIF 


И 
LINE(K ) =KODE( 1 
1 CONTINUE 
NPLOT = NPLOT + 1 
WRITE(6,1) 
0 CONTINUE 


FORMAT (1X) 

FORMAT (1X,217) 

FORMAT (1X,14,3F9.4,F10.0) 
FORMAT (1X,F6.4,4F8.4,90A1) 


END 


7439/7 


NUMDS: NPLOT: 
2 So 


»CPY(2,L),CPY(3,L),LINE 


CPU: 
CPY(2,L).CPY(3,L),LINE 


ОЩЕСЕ ВЕ 


РКОСКАМ СЕТХ 


THIS PROGRAM READS COMBINED (MULTIPLE DATA SET) FILES 
OF XYZ AND Q PLOTTING DATA AND SEPARATES THEM INTO 
INDIVIDUAL FILES FOR PLOTSD SUBMITTAL. 


“ 
° 
$ 
° 
° 
VARIABLES: ° 
CFXNAME: NAME OF COMBINED XYZ FILE 4 
CFONAME: NAME OF COMBINED Q FILE . 
FXNAME: NAME OF INDIVIDUAL X FILES, CHAR VARIABLE ° 
FONAME: NAME OF INDIVIDUAL Q FILES, CHAR VARIABLE 2 
DSSKIP: NUMBER OF DATA SETS TO SKIP WHEN à 
READING COMBINED FILE з 

DSREAD: NUMBER OF DATA SETS TO READ з 
INTVL: INTERVAL BETWEEN DATA SETS SENT TO FILE ° 

a 

a 


O. 8e 8 а а а са е а с е е е е е е га С) 


INTEGER DSSKIP,DSREAD,COUNT 

CHARACTER FXNAME + 12, FONAME « 12, ,CFXNAME + 10, CFONAME + 10 

DIMENSION X(250,100),2(250,100) 

DIMENSION 01(250,108),02(250,100),03(250,100),04(258,100) 
С 
eee USER INPUT: е«еееееегевеегоеоееееееееееееееееевеееегогееееевеееее 
С 

DATA XSOGALE, YSCALE/1.,1./ 

DATA XPOSIT,YPOSIT/-1.,-1./ 

DATA DSSKIP/10/ , DSREAD/3/ , INTVL/1/ 

СЕХМАМЕ = '38x.iris' 

CFQNAME = '38q.iris' 
С 
0.0.0 99. 1.9999 099990 099 1 999090 1.99 09 99 999 9 90 99 99 99 19 119 909965 
C 

OPEN(UNIT=91,FILE=CFXNAME , STATUS="OLD” , FORM='BINARY') 

OPEN(UNIT=92,FILE=CFONAME , STATUS="OLD” ,FORM="'BINARY? ) 


FXNAME = 'x000. iris” 
FONAME = 'q000. iris” 
COUNT = 0 


IF (DSSKIP.GT.@) THEN 
DO 1@ ISKIP = 1,DSSKIP 
READ(91) IMAX,KMAX 
КЕАО(91) ((Х(І,К), 1=1, ІМАХ),К=1,КМАХ), 
+ ((ZCI,K),Iz1, IMAX) ,Kz1,KMAX) 
READ(92) IMAX,KMAX 
READ(92) A,B,C,D 
READ(92) ((Q1(1.K),1=1, IMAX),K=1,KMAX), 
((02(1,К), 1=1, | МАХ) ,К<1,КМАХ), 
((03(1,К), 11, 1МАХ ) ,K=1 . KMAX ) , 
((Q4(1,K),1=1, IMAX) ,K=1,KMAX) 


+++ 


10 CONTINUE 
ENDIF 


DO 20 IREAD = 1,DSREAD 
COUNT = COUNT + 1 
WRITE (FXNAME(2:4),100) IREAD 
WRITE (FONAME(2:4),100) IREAD 
READ(91) IMAX,KMAX 
READ(91) ((Х(1.К).1=1,1МАХ).К=1 , КМАХ), 
> ((Z(1,K). I=1, IMAX) ,K=1,KMAX) 
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12 
11 


190 


+++ 


+++ 


IF (COUNT .EQ. INTVL) THEN 
OPEN(UNIT=1,FILE=FXNAME, FORM='BI]NARY ” ) 
DO 11 Iz1,1MAX 

DO 12 K=1,KMAX 
X( 1 ,K)=X( 1 ,K) sXSCALE+XPOS 1T 
Z(1,K)=Z(1,K)*YSCALE+YPOSIT 
CONT I NUE 
CONTINUE 
WRITE( 1 ) IMAX,KMAX 
WRITE( 1 ) ((Х(Т,К), 1=1, [МАХ) ‚К=1,КМАХ), 
((2(1.К), 11, МАХ) ‚К=ж1,КМАХ ) 


CLOSE( 1 ) 
ENDIF 
READ(92) IMAX,KMAX 
READ(92) A,B,C,D 
READ(92) ((Q1(1,K),I=1,]MAX),K=1,KMAX), 
((02(1,К),1=1, ТМАХ ) „К=1 ,КМАХ ), 
((03(1,К), 11, ТМАХ ) ,K=1 . KMAX ) , 
((04(1,К),1=1, ТМАХ ) „К=1 ,КМАХ ) 
IF (COUNT.EQ. INTVL) THEN 
OPEN(UNIT=2,FILE=FQNAME , FORM='BINARY ' ) 
WRITE(2) IMAX,KMAX 
WRITE(2) A,B.C.D 
WRITE(2) ((Q1(1,K),I=1,IMAX),K=1,KMAX), 
((Q2(1,K), 1=1, IMAX) ,K=1,KMAX), 
((Q3(1,K),1=1, IMAX) ,K=1,KMAX), 
((04(1,К), 11, ТМАХ ) ,K=1,KMAX) 
CLOSE( 2 ) 
COUNT = 0 
ENDIF 
CONTINUE 
CLOSE(91) 
CLOSE(92) 
FORMAT( 13.3) 
STOP 
END 


pa GERCTPF 


PROGRAM GETCP 


$ 

° THIS PROGRAM ACCESSES COMBINED CP OR CF FILES ON THE ° 
° IRIS ACCOUNT, SEPARATES THEM INTO INDIVIDUAL PLOTTING ә 
ゅ FILES AND SCALES THEM TO MEET SPECIFIED REQUIREMENTS. > 
° A DUMMY Q FIL* :S ALSO GENERATED (Q. IRIS). ° 
Ф 3 
в VARIABLES: 5 
ອ Ф 
° CPNAME: COMBINED FILE NAME e 
в ЕМАМЕ: OUTPUT FILE NAME, CHARACTER VARIABLE » 
° DSSKIP: NUMBER OF DATASETS IN THE COMBINED FILE е 
> TO SKIP ° 
в DSFILE: NUMBER OF DATASETS IN THE COMBINED FILE a 
° TO SEPARATE AND FILE ə 
ゅ XSCALE: AXIS SCALING FACTOR > 
+ YSCALE: CP/CF SCALING FACTOR (RELATIVE MAGNITUDE) ә 
° XPOSIT: X POSITION FACTOR (SCREEN LOCATION) > 
ゅ YPOSIT: Y POSITION FACTOR (SCREEN LOCATION) в 
~ # 
Ф = 
C 


INTEGER DSSKIP,DSFILE 
CHARACTER ҒМАМЕ • 1 1, ОМАМЕ •6 , СЕМАМЕ • 10 
DIMENSION X(3,75).2(3,75),XS(3.75).25(3,75) 
DIMENSION Q1(3,75),02(3,75),03(3,75),04(3,75) 
ОПАТА 01,02,05,04/900•1./ 


C 
ese USER INPUTS: зеоввваовозововово овоа ово 
し 
DATA DSSKIP/256/ , DSF ILE/44/ 
DATA XSCALE/.5/.YSCALE/-.05/ 
DATA XPOSIT/-1.0/,YPOSIT/.5/ 
СЕМАМЕ = 'snkrO5p' 
FNAME = ‘'cp0@@.iris’ 
C 
s6252530999909950995090995950995990599999565999990995959599995059529$$909 
C 
OPEN(UNIT«90,FILEZCFNAME, STATUS&'OLD' , FORME'BINARY ') 
C 
ОНАМЕ ж 'q.iris' 
С 
eee SKIP DSSKIP FILES 
C 
IF (DSSKIP.GT.0) THEN 
DO 10 ISKIP = 1,DSSKIP 
КЕАО (90) ІМАХ,КМАХ 
READ( 99 ) ((Х(1,К) , 1=1, IMAX) , Kx1, KMAX) , 
+ ((2(1,К), 1=1,1МАХ),К=1 ‚КМАХ) 
10 CONTINUE 
ENDIF 
C 


sss SEPARATE AND SCALE DSFILE FILES 


DO 40 IFILE = 1,DSFILE 
WRITE(FNAME(3:5),100) IFILE 
OPEN(UNIT=1, FILE=FNAME, FORM=’BINARY' ) 
READ(9@) IMAX,KMAX 
READ(90) ((X(I,K),I=1,IMAX),K=1,KMAX), 


((Z(1,K),1=1, IMAX) ,K=1,KMAX) 


29 
39 


199 


DO 39 I = 1,IMAX 
DO 29 K = 1,KMAX 


XS(I,K)=X(I,K)*XSCALE + XPOSIT 
ZS(1,K)=Z(1,K)»YSCALE + YPOSIT 


CONTINUE 
CONTINUE 
WRITE(1) IMAX,KMAX 


WRITE(1) ((XS(1,K), I=1, IMAX) ,Ка1 , КМАХ 
+ ((25(1,К), 1=1, ТМАХ ) „К=1 ,КМАХ 


CLOSE( 1 ) 
CONT INUE 


GENERATE DUMMY Q FILE 


OPEN(UNIT=1,FILE=QNAME , FORM='BINARY' ) 
WRITE(1) МАХ, КМАХ 

WRITE( 1 ) IMAX, IMAX, IMAX, 1МАХ 

WRITE( 1 ) A TEM MEN ME 


+ ((Q2(1,K),1=1, IMAX) .K=1,KMAX), 

+ ((Q3(1,K),1=1, IMAX),K=1,KMAX), 

+ ((04(1,К),1=1, ТМАХ ) ,Ke1,KMAX) 
CLOSE( 1 ) 


FORMAT(13.3) 


STOP 
END 
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8. GETSIN.F 


PROGRAM GETSIN 


ゆり りり おち りり りり わり ゆり わり ゆり りり りり あり りあ お あお お お あり りあ ああ ちち ちり お りり りり $$ 


の 
° THIS PROGRAM READS A COMBINED FILE OF AOA POINTER 
в ПАТА (СЕМЕКАТЕО WITH PLOT3D FUNCTION 0 ) AND SEPARATES 
ә IT INTO INDIVIDUAL PLOTTING FILES. 
4 VARIABLES: 
P CFNAME: NAME OF COMBINED FILE 
“ ҒМАМЕ: NAME OF INDIVIDUAL FILES, CHAR VARIABLE 
- DSSKIP: NUMBER OF DATA SETS TO SKIP WHEN 
* READING COMBINED FILE 
+ DSFILE: NUMBER OF DATA SETS TO SEND TO INDIVIDUAL 
ж FILES 
° 
ооо оо оо ооо ооо оо ооо ово о ооо ооо оо оо ооо ооо оо ооо ооо оо ооо ооо 
С 
INTEGER OSSKIP,DSFILE 
CHARACTER FNAME + 12 ,CFNAME +10 
DIMENSION X(3,294),2(3,294) 
C 
eee USER INPUT: eeeee ee の の の ひ の の の の の の の の の の の の の ゆ ゆ の の ゆい の の の の の の の ゆ ゆ の の の ゆ ゆ の の の の の ゆ の 
Ф 
DATA DSSKIP/250/,DSFILE/44/ 
СЕМАМЕ = 'впкгд55' 
С 
ооо ооо оо ооо соо ооо осо ооо ооо ооо ооо о ооо ооо оо ооо оо ооооо 
C 
OPEN(UNIT=90,FILE=CFNAME , STATUS=*OLD” , FORM="BINARY” ) 
C 
FNAME = 'sin@@@.iris’ 
С 
IF (DSSKIP.GT.@) THEN 
DO 10 ISKIP = 1,DSSKIP 
READ(90) | МАХ, КМАХ 
READ( 99 ) ((Х(1,Ј), 121, ТМАХ ) ,Ј=1 ,КМАХ ), 
+ ((2(1,Ј) , 11, ]ЈМАХ) , Ј=1 ,„КМАХ ) 
19 CONTINUE 
ENDIF 
Ф 
DO 28 IFILE = 1,DSFILE 
WRITE (FNAME(4:6),100) IFILE 
OPEN(UNI T=1,FILE=FNAME , FORM='BINARY') 
READ(9@) ІМАХ,КМАХ 
READ( 99 ) ((Х(1,Ј), 1=1, ТМАХ) , J=1 .KMAX ) 。 
+ ((2(1,7Ј), 1=1, ІМАХ), Ј=1, КМАХ) 
WRITE( 1 ) IMAX , KMAX 
WRITE(1) ((Х(1,9), 1=1, ТМАХ ) , ЈЕ1 ,КМАХ ),, 
+ ((Z(1,J),I=1,IMAX),J=1,KMAX) 
CLOSE( 1 ) 
20 CONTINUE 
С 
CLOSE(90) 
100 FORMAT(13.3) 
STOP 
END 
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C. PROCEDURAL DOCUNIENTATION 


After proper initialization (user-inputs) as described above, all codes or JCL s which 


produce or transfer Plot3D files are subnutted to the Crav utilizing standard CSUB 


commands. Once the data is resident on the ERIS and tlie user is logged on and in the 


proper directory level. the following procedures may be utilized for graphics gancration. 


l. 


ry 


Separate the combined file into discrete data sets. 


Edit “user inputs” section of GETA.F (or other files as ipnropriatey io (ei eee 
me (ultimate location on the displav), and identify the conibined feto асс 
those «ata sets to De ses ще 


vi getx.f 

Сошруе the editted program: 
(77 getx.f 

un the compiled version bv entering the filename (2.0ut by defauit). 
a.out 


Discrete data sets will appear on the account with names corresponding to those 
in the Plot3D initialization files. Check by using the Is command. 


(зе Plot3D to generate ora iien 


Use the mv commands to give the proper initialization Ше Шо Са 
“plot3dini.com”. (fhe [RIS will not save multiple versions of files with the same 
filename. Instead, older versions of the file will be deleted 15 Order to NE 
advertent deletion of initialization files, multiple mv conimands may be required.) 


Inv Xini.com plot3dini.com 


If processing flow field files (X and Q), edit the initialization file to identify desired 
function and number of contours. This involves changing only two lines at the top 
Sl one Nie. 


vi plot3dini.com 

Create era ШЕЕ 
plot3x 

Rename the resultant combined graphics file. 
mv q001.gra somefilnm.gra 


After creation of as manv .gra files as required froin the separated files, clean up 
the account using the following wildcards. 


rm q*.1ris 


rm X*. iris 


48 


ຈ 
Js 


[nteractive viewing. 


Run the Graphies Animation System software. 


048 


> 


Input .gra [iles to GAS bv selecting on the main menu: 


ARCGRAPH file input 


On the submenu. select: 


load entire file 


[n response to prompts. enter the sequential object number. aud 


( someflnm.gra ). 
colormap. 


1 


yer color plots, press БЕРСКА when prompt 


Me the data by selecting on the main menu: 


view data 


Ltilize the mouse to manipulate the display. 
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Í 
ed 


APPENDIX B. SANKAR NAVIE STOKES SOLVER 


(JCL) 


/EOF 
PROCRAM MAIN 


SANKAR NAVIER STOKES CODE: SOLVES TWO DIMENSIONAL FLOWS PAST 
ARBITRARY GEOMETRIES USING ADI PROCEDURE. THIS VERSION SAVES 
MULTIPLE PLOT3D DATA SETS IN SINGLE FILES. 


Ф 
е 
Ф 
е 
а 
VARIABLES: * 
TITLE: (60 CHARACTERS MAX.) з 
IMAX: NUMBER OF X COORD LOCATIONS 3 
KMAX: NUMBER OF Y COORD LOCATIONS * 
ITEL: GRID POINT AT LOWER TRAILING EDGE з 
ITEU: GRID POINT AT UPPER TRAILING EDGE e 
ОТ: SIZE OF TIME STEP 3 
WW: EXPLICIT ARTIFICIAL VISCOSITY TERM 3 
ALFA: MEAN ANGLE OF ATTACK (DEGREES) 3 
ALFA1: AMPLITUDE OF OSCILLATION (DEGREES) * 
ALFAI: ANGLE BELOW WHICH MODIFIED TURBULENCE MODEL USED TO * 
COMPUTE EDDY VISCOSITY * 

REDFRE: REDUCED FREQUENCY d 
AMINF: FREE STREAM MACH NUMBER > 
REYREF: REYNOLDS NUMBER (MILLIONS; NEG. = INVISCID FLOW) . 
DNMIN: DISTANCE OF FIRST POINT OFF OF WALL a 
TSTART: TIME FOR CALCULATIONS TO START ON PRESENT RUN ຈ 
FORMAT: OLDSLN FORMAT ;3.0=PLOT3D FILES,-3.0=Q MATRICES з 
RSTRT: TRUE = STORED DATA READ TO CONTINUE ITERATION 9 
PITCH: TRUE = AIRFOL AOA OSCILLATES a 
PLUNGE: TRUE = AIRFOIL OSCILLATES VERTICALLY з 
FNU: NUMBER OF UPPER SURFACE DEFINITION POINTS ° 
FNL: NUMBER OF LOWER SURFACE DEFINITION POINTS * 
FSYM: SYMMETRY FLAG (1 = SYMMETRIC) 3 
X: AIRFOIL GEOMETRY DEFINITION POINTS ә 
Y: AIRFOIL GEOMETRY DEFINITION POINTS a 
CSTP: NUMBER OF STEPS COMPLETED 3 
CPLT: NUMBER OF DYNAMIC PLOTS COMPLETED в 
NSTP: NUMBER OF TIME STEPS TO BE COMPLETED ON PRESENT RUN @ 
PSTP: NUMBER OF TIME STEPS BETWEEN PLOT DATA OUTPUT э 
$ 

@ 

d 

в 

Ф 

a 

Ф 

@ 

a 

a 

@ 

a 


TAPE DEFINITIONS: 
TAPE 05: INPUT DATA 
TAPE 06: OUTPUT DATA 
TAPE 07: FLOW FIELD INPUT DATA FOR RESTARTS 
TAPE 29: PLOT3D XYZ FILE STORAGE 
TAPE 21: PLOT3D Q FILE STORAGE 
TAPE 50: PLOT3D CP XYZ FILE STORAGE 
TAPE 60: PLOT3D CF XYZ FILE STORAGE 
TAPE 90: PREVIOUS RUN X FILE STORAGE 
TAPE 91: PREVIOUS RUN Q FILE STORAGE 
TAPE 92: PREVIOUS RUN CP FILE STORAGE 
TAPE 93: PREVIOUS RUN CF FILE STORAGE 


тгееесесеееееоееееееееоееееееоеееееееееееееее;ееееееееееееееевееееефеееееееевееееее 
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INTEGER PSTP,CSTP,CPLT 
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Cees 


Cees 


ດດດ 


Cece 


COMMON/ SURF /PSUR ( 161) 

COMMON/FIX/OMEGA 

COMMON/MUTUR/CMU( 161,41) 

COMMON/ SK INCF /CF ( 161) 

COMMON/GRID1/X(161,41),2(161,41) 

COMMON /PAR/GAMMA , REYREF , ALFA, ALFAT,REDFRE,AMINF.,ALFAI 
COMMON/DGRID/DT, IMAX,KMAX,ITEL,ITEU 

COMMON/GRID/YACOB( 161,41) 

COMMON /DAMP /WW , WW I 

COMMON/FLOW/Q1( 161,41) ,Q2(161,41),Q3(161,41),Q4(161, 41) 
DIMENSION TITLE(15) ,TYTLE(15) . ALPHA( 96 ) ,CPTH( 97 , 96 ) . XTH( 97 ) 
DIMENSION XR(161,49),ZR(161,49),01R(161,41),02R(161,41) 
DIMENSION Q3R(161,41),Q4R(161,41) 

COMMON/ LOGIC/RSTRT,PITCH,PLUNGE 

LOGICAL RSTRT,PITCH,PLUNGE 


READ INPUT DATA 


READ N 

READ (5,1) TITLE 

READ с 

READ (5,2) IMAX,KMAX,DT,WW,ALFA,ALFA1,ALFAI,REDFRE,AMINF 

READ (5,5 

READ (5,3) ITEL, ITEU, REYREF, ONMIN, TSTART, FORMAT, RSTRT,PITCH,PLUNGE 
READ (5,5 

READ (5,4) CSTP,CPLT.NSTP,PSTP 


INITIALIZE 


GAMMA = 1.4 

РТ = АТАМ(1. )•4. 

SET ALFAD FOR STEADY STATE PLOT3D OUTPUT 

ALFAD = ALFA 

FORCE DT TO BE EQUAL TO UNITY FOR STEADY FLOW PROBLEMS 
THIS INVOKES THE RELAXATION MODE 

IF(REDFRE.LE.0.001) DT = 1.0 

REYREF = REYREF * 1.E+@6 

NITN = CSTP + NSTP 

NPLOTS = CPLT 

CPLT = CPLT + 1 

CSTP = CSTP + 1 

DENSITY NORMALISED WITH RESPECT TO ROINF 

VELOCITIES NORMALISED WITH RESPECT TO AINF 

TOTAL ENERGY NORMALISED WITH RESPECT TO (ROINF*AINF*AINF) 
ТОТЕ№=АМ 1 НЕ АМ 1МЕ*0. 5+1. / (САММА_* (GAMMA-1.) ) 

ALFA = ALFA = ATAN(1.) / 45. 

ALMEAN = ALFA 

ALFAI = ALFAI * ATAN(1.) / 45. 


ALFA1 = ALFA1 * ATAN(1.) / 45. 

ALFAS = ALMEAN — ALFA1 

WRITE(6)' PLOT ITN TIME AOA CH 60 СМ” 
UINF = AMINF 

VINF = 0.0 


DO 7 I=1,IMAX 

DO 7 K=1,KMAX 
Q1(1,K)=1. 
Q2(1,K)=UINF 
03(1,K)=VINF 
04(1,K)=TOTEN 

CONT INUE 


INPUT STORED FLOW FIELD DATA 


IF(RSTRT )THEN 
IF (FORMAT.LT.@.@) THEN 
READ(7) TIME,01,02,03.04 
ELSE 
READ(7) ІМАХ,КМАХ 


5] 


READ(7) AMINF ,ALFAD,REYREF , TIME 

READ(7) ((Q1(I,J), I=1,IMAX),J=1,KMAX), 
((О2(1,Ј), I=1,IMAX),J=1,KMAX), 
((Q3(1,J), I=1,IMAX),J=1,KMAX), 
((Q4(1,J), I=1,IMAX),J=1,KMAX) 


+++ 


ENDIF 
IF(NPLOTS.GT.@)THEN 
DO 9 Ki=1,NPLOTS 
КЕАО (98) IMAXR,KMAXR 
КЕАО (99) ((XR(1,J), I=1,IMAXR),J=1,KMAXR), 
+ ((ZR( 1 ,J) , I=1,IMAXR),J=1,KMAXR) 
WRITE(20) IMAXR . KMAXR 
WRITE(20) ((XR(I,J), 1=1 , IMAXR ) , J=1 . KMAXR ) , 
+ ((ZR(I,J), I=1,IMAXR) ,J=1,KMAXR) 
9 CONTINUE 
DO 1 の K2=1 .NPLOTS 
READ(91 )  IMAXR , KMAXR 
READ(91) AMINFR,ALFAOR,REYREFR, TIMER 
READ(91) ((Q1R(1,J), „ ТМАХК ) , Ј=1 ,КМАХЕ ), 
,IMAXR) , J=1 , KMAXR ) , 


[21 

+ ((02R(1,J), I=1 
гр ((Q3R(1,J), I=1,IMAXR),J=1,KMAXR), 
S ((Q4R(I,J), I21, IMAXR) , J21, KMAXR) 

WRITE(21) IMAXR,KMAXR 

WRITE(21) AMINFR,ALFADR,REYREFR, TIMER 

WRITE(21) ((QIR(I,J), 1=1. IMAXR ) , J=1 , KMAXR ) , 
+ ((Q2R(I,J), I=1,IMAXR),J=1,KMAXR), 
+ ((Q3R(I,J), I=1,IMAXR),J=1,KMAXR), 
+ ((Q4R(I,J), I=1,IMAXR),J=1,KMAXR) 

10 CONT INUE 


DO 11 K3=1,NPLOTS 
READ(92)  IMAXR,KMAXR 
READ(92) ((XR(1,J), I=1,IMAXR),J=1,KMAXR), 
- ((ZR( 1 .J) , 1=1 , IMAXR ) , J=1 , KMAXR ) 
WRITE(S59 ) IMAXR , KMAXR 
WRITE(58) ((XR(1,J), 1!=1 , IMAXR ) , J= 1 .KMAXR ) , 
+ ((ZR(I,J), I=1,IMAXR),J=1,KMAXR ) 
11 CONTINUE 
DO 12 K4=1.NPLOTS 
READ(93)  IMAXR,KMAXR 
READ(93)  ((XR( 1 ,J) , 1=1 , IMAXR ) , J=1 , KMAXR ) , 
+ ((ZR(I,J), 1=1, МАХВ) , J=1,KMAXR) 
WRITE(60) IMAXR,KMAXR 
WRITE(60) ((XR(I,J), I€1, IMAXR) , J21, KMAXR), 
+ ((ZR(1,J), I=1,IMAXR),J=1,KMAXR) 
12 CONTINUE 
ENDIF 
ENDIF 
IF(TSTART.GE.@.) TIME = TSTART 
IF(.NOT.(RSTRT)) TIME = 0. 


C 

Се. GENERATE COMPUTATIONAL GRID, DEFINE SURFACE COORD © © AOA, 
C ROTATE TO INITIAL AOA AND COMPUTE METRICS 

C 


CALL AIRFOL(IMAX,KMAX,ITEL,ITEU) 
IF (REYREF.GT.@) CALL CLUSTR(DNMIN) 
ТЕОСЕ = X(ITEL + 48,1) 
DO 15 Tæ 1,97 
XTH(1) = X(I+ITEL-1,1)-TEDGE 
ГЭ CONTINUE 
CALL ROTGRID(X,Z, IMAX,KMAX , ALFAS) 
CALL METRIC 
Ç 
Csss ITERATIVE LOOP 
C 


00 1000 ІТМе1 ,М5ТР 
TIME = TIME + DT 


Cees ROTATE GRID TO NEW ANGLE OR SET TO CORRECT ANGLE FOR RESTARTS 
C 
IF (PITCH) THEN 
OMEGA = 2. • REDFRE *AMINFeSIN(REDFRE * 2. * TIME * AMINF) 


1 ФА| РА! 
ALOLD = ALMEAN - ALFAI = COS(2. = REDFRE = AMINF з 
1 (TIME - DT)) 


ALFA = ALMEAN - ALFA1 + COS(REDFRE • 2. • ТІМЕ e AMINF) 
ALFAD = ALFA • 45. / АТАМ(1.) 
DALFA = ALFA — ALOLD 
IF(RSTRT.AND. ITN. EQ. 1) DALFA = ALFA — 2.s+ALFAS 
CALL ROTGRID(X,Z, IMAX,KMAX,DALFA) 
CALL METRIC 
END IF 
IF (PLUNGE) THEN 
OMEGA = 2. * REDFRE * AMINF 
ALFA = ALMEAN 


END IF 
C 
Cese COMPUTE THE SOLUTION BY ADI TECHNIQUE 
С 
CALL SLPS(ITN,OMEGA,DALFA) 
し 
Севе APPLY BOUNDARY CONDITIONS 
6 
CALL WALLBC 
C 
Cees GENERATE PLOT3D FILES 
С 
IF(CSTP.EQ.(CPLTePSTP)) THEN 
WRITE(20) IMAX, KMAX 
WRITE(20) ((X(1,J), I=1,IMAX), J=1,KMAX), 
1 ((Z(1,J), J=1,IMAX), J=1,KMAX) 
WRITE(21) IMAX, KMAX 
WRITE(21) AMINF, ALFAD, REYREF, TIME 
WRITE(21) ((Q1(I,J), Iz1,IMAX), Jz1,KMAX), 
1 ((О2(1.Ј), I=1,I]IMAX), J=1,KMAX), 
2 ((ОЗ(1,Ј), Ге1 ,ТМАХ), J=1,KMAX), 
3 ((О4(1,Ј), Т=1,]МАХ), Ј=1,КМАХ) 
С 
Cees GENERATE PERFORMANCE COEFFICIENTS 
C 


CALL LOAD(PSUR ,CL , CD ,CM, ALFAS ) 
AOA = ALFA*189./PI 
WRITE(6.6) CPLT,CSTP,TIME,AOA,CL,CD,CM 
CALL CPPLOT(ITEL, ITEU, AMINF,X(1,1),2(1, 1) , PSUR) 
CPLT = CPLT+ 1 

END IF 

CSTP = CSTP + 1 

1000 CONTINUE 


FORMAT (1X,15A4) 

FORMAT (1X,217,7F8.4) 

FORMAT (1X,217,4F8.5,3L7) 

FORMAT (1X,417) 

FORMAT (1X) 

FORMAT (1X,13,16,F8.3,F9.5,3F8.4) 


END 


OOO OD の On + Gi Ñ — O 


SUBROUTINE AMAT1(K,IMX1,XIX,XIZ,XIT) 
COMMON/FLOW/Q1(161,41),02(161,41),03(161,41),04(161,41) 
COMMON/PERTR/DQ1(161,41),D02(161,41),D03(161,41),D04(161,41) 
COMMON /AM/A( 4.4, 161) 

COMMON /PAR/GAMMA ,REYREF ,ALFA,ALFA1,REDFRE, AMINF ,ALFAI 


EE 


Севе 


1999 


C 


DIMENSION XIT(161,41),XIX(161,41),XIZ(161,41) 
REAL K0,K1,K2 


AMAT1 COMPUTES THE COEFFICIENT MATRIX DE/DQ DURING XI SWEEP 


GM1 = GAMMA - 1. 

00 1000 I = 2 , IMX1 

ке = XIT(I,K) 

K1 = XIX(I,K) 

K2 = XIZ(I,K) 

U z Q2(I,K) / Q1(I,K) 

W «, ОЗ(І,К) / ОТ(І,К) 

EBYR = Q4(1,K) / Q1(1,K) 

PHI2 = 0.5 e GMI e (Us U+ W» W) 

THETA = Ki * U + K2 * W 

A(1,1,1) = K@ 

А(1,2,1) = КЇ 

T = K2 

A(1,4,1 = Q 

A(2,1,1) = K1 * PHI2 — U * THETA 

A(2,2,1) = КӨ + THETA — K1 e (CM1 - 1.) e U 
A(2,3,1) = K2 * U — СМ1 е КТ ей 

А(2,4,1) = K1 * GM1 

A(3,1,1) = K2 * PHI2 — W e THETA 

A(3,2,1) = K1 * W — K2 * GM1 e U 

А(3,3,1) = KQ + THETA — K2 e (GM1 — 1.) * W 
А(3,4,1) = K2 * GM1 

A(4,1,1) = THETA e (2. e PHI2 — GAMMA = EBYR) 
А(4,2,1) = K1 = (GAMMA » EBYR - PHI2) — GM1 s U èe THETA 
А(4,3,1) = K2 * (GAMMA » EBYR - PHI2) — GM1 e W «= THETA 
А(4,4,1) = К@ + GAMMA * THETA 

CONTINUE 

RETURN 

END 
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C 


Сезе 


SUBROUTINE AMAT2(1,KMX1,ZETAX, ZETAZ,ZETAT) 

COMMON /FLOW/Q1(161,41),Q2(161,41),Q3(161,41),04(161, 41) 
COMMON/PERTR/0Q1(161,41),002(161,41),003(161,41),004(161,41) 
COMMON/AM/A(4,4,161) 
COMMON/PAR/GAMMA , REYREF , ALFA, ALFA1, REDFRE, AMINF,ALFAI 
DIMENSION ZETAX(161,41),ZETAZ(161,41) , ZETAT(161,41) 

REAL KO,K1,K2 


AMAT2 COMPUTES THE COEFFICIENT MATRIX DF/OQ DURING ETA SWEEP 


GM1 = GAMMA - 1. 
ОО 1000 К = 2 , KMX1 

ке = ZETAT(I,K) 

K1 = ZETAX(I,K) 

K2 = ZETAZ(I,K) 

U ж 02(І,К) / 01(І,К) 

W z Q3(I,K) / Q1(I,K) 

EBYR = Q4(I,K) / Q1(1,K) 

PHI2 = Q.5 * GM1 * (U * U + W е W) 
THETA = K1 e U + K2 * W 

A(1,1,K) = K@ 

А(1,2,К) = K1 

A(1,3,K) = K2 

A(1.4.K) = ө 

A(2,1,K) = K1 * PHI2 -U * THETA 
A(2,2,K) = KO + ТНЕТА - К1» (СМ1-1.) ей 
A(2,3,K) = K2 * U — GMI е КТ е И 
A(2.4.K) =K1 * GM1 

А(3,1,К) -К2 е РНІ2- № « THETA 
A(3,2,K) = K1 е W- K2 © GM1 e U 
A(3,3,K) = KQ + THETA —K2 * (GM1-1.) * W 


K2 * GM1 

THETA • (2. «e PHI2 — GAMMA « EBYR) 

K1 e (GAMMA « EBYR — PHI2) — СМ! • Џ • ТНЕТА 
К2 + (GAMMA e EBYR — PHI2) — GM1 e W e THETA 
KO + GAMMA e THETA 


M3, ,4,K) 
A(4,1,K) 
າ 
A(4,3,K 
A(4,4,K) 

1000 CONTINUE 
RETURN 
END 


C 
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C 
SUBROUTINE SLPS(ITN,OMEGA,DALFA) 
COMMON /FLOW/Q1(161,41),02(161,41),03(161,41),04(161,41) 
COMMON/PERTR/DOQ1( 161,41) ,002(161,41) ,003(161,41) ,004(161,41) 
COMMON/AM/A(4,4,161) 
COMMON/TRID/DD(4,4,161,41),MM(4,4,161,41) ,EE(4,4,161,41) 
1,GG(4,161,41) 
COMMON /PAR/GAMMA , REYREF , ALFA, ALFAT , REDFRE,AMINF , ALFAI 
COMMON/DGRID/DT, IMAX,KMAX, ITEL, ITEU 
COMMON/GRID/YACOB(161,41) 
COMMON/DAMP /WW , WW 1 
COMMON/MTRIX/ X1X(161,41),XIZ(161,41).ZETAX(161,41), 
*ZETAZ(161,41) 
1 ,XIT(161,41),ZETAT(161,41) 
REAL MM 
DIMENSION DELTAT(161,41) 
C 
Cees SUBROUTINE SLPS DOES THE BULK OF THE WORK FOR THE AL! ALGORITHM. 
Ceee IT CALLS FLUX AND COMPUTES RIGHT HAND SIDE DURING THE TWO SWEEPS, 
Case ASSEMBLES THE COEFFICIENT MARICES, ADDS IMPLICIT AND EXPLICIT 
Ceee DISSIPATION AND CALLS THE TRIDIAGONAL SOLVER TO OBTAIN THE FINAL 
Ceee SOLUTION. 
C!!!!iSET VALUE OF IMPLICIT DAMPING COEFFICIENT 
WWI = 20.0 +. WW 
IM1 = IMAX 一 1 
IM2 = | МАХ 一 2 
KM1 = KMAX 一 1 
KM2 = KMAX — 2 
IF(ITN.EQ.1) THEN 
IF(REDFRE.LT.0.001) THEN 
DO 777 K = 2 , KMAX — 1 
DO 777 1 = 2 , IMAX - 1 
ОЕ(ТАТ(І,К) - 9.5 / (1. + SORT(ABS(YACOB(I,K)))) 
777 CONTINUE 
ELSE 
ОО 778 K = 2 , КМАХ - 1 
DO 778 I = 2 , МАХ - 1 
778 DELTAT(1,K) = 1.0 
END IF 
END IF 
C 
Ces» THE DISSIPATION TERMS ARE CONSTRUCTED AND STORED IN THE ARRAYS 001, 
Csse DQ2,DQ3 AND 004. 
し 
CALL DISSIP 
C 


Cese THE RIGHT HAND SIDE AT KNOWN TIME LEVEL IS NOW COMPUTED AND ADDED 
CALL RESI (OMEGA) 
E 
Ceee JF VISCOUS FLOW IS COMPUTED THE VISCOUS TERMS ARE ADDED TO DQ1 ETC. HERE. 
С 
IF(REYREF.GT.@.) CALL STRESS(ITN,DALFA) 
Севе 1-5МЕЕР. 
し 
DTH = DT e 0.5 
DTW = DT • WWI 
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O O О 


OOO 


990 


991 


17 
13 


DO 3 K 


2 , КМ 


CALL AMAT1(K,IMAX-1.XIX,XIZ.XIT) 


DO 4 [1 
DO 4 [2 
DO Š 1 


ЕЕ(І1.12,1-1,К) 
00(11.,12,1-1,К) 


CONT INUE 
CONTINUE 


1 , 4 
1 , 4 
2. ІМАХ - 1 


А СЕТ. 12 репу < DTH ПЕТАК) 
– А(11,12,1-1) • ОТН • ОЕ(ТАТ(1„К) 


IMPLICIT DAMPING ADDED HERE. 


DO 6 11 
DO 6 1 
DT1 


DD(11,11,I-1,K) 
EE(11.11,1-1,K) 
‚11, [-1,К) 


мм(11 
CONT I NUE 
DO 990 I 


66(1.1-1,К) 


GG(2,I-1,K 
GG(3,I-1,K 


| 


СО К) 


CONTINUE 
CONTINUE 


1,4 
2 , [IMAX =- 1 

DTW / YACOB(I,K) + DELTAT(I,K) 
DD(I1.11,I-1,K) - DT1 * YACOB( I-1,K) 
ЕЕ(11,11,1-1,К) - ОТ! г YACOB(I*1,K) 
1. + 2. e DIW e DELTAT(I,K) 


2 , IMAX — 1 
DO1(I,K) * DELTAT(I,K) 
DQ2(I,K) e DELTAT(I,K) 
DO3(I,K) e DELTAT(I,K) 
004( ! ,“) * DELTAT(I,K) 


PERFORM BLOCK TRIDIAGONAL MATRIX INVERSION FOR THE ENTIRE PLANE 


CALL MATRX1(IMAX,KMAX ) 


DO 991 K 
00 991 1 
DQ1(I.K) 
DO2(1,K) 
DO3(I.K) 
DO4(I.K) 
CONT INUE 


2 , KMAX - 1 
2 , IM1 

GG(1,I-1,K) 
GG(2,I-1,K) 
GG(3,I-1,K) 
GG(4,I-1,K) 


K-SWEEP BEGINS HERE. 


DO 13 ! = 2 , IM1 
CALL AMAT2(I,KMAX-1,ZETAX,ZETAZ,ZETAT) 
ВО он 1 = 1, 4 

DO 15 12 72 

DO 15 K 2 , KMAX - 1 


EE(11.12,1,K-1) 
рр(11,12,1,К-1) 


CONTINUE 


A(I1,12,K+1)*OTH e DELTAT(I,K) 
—-A(I1,I2,K-1)*DTH * DELTAT(I,K) 


SECOND ORDER DAMPING ADDED HERE. 


DO 16 11 
DO 16 K 
DT1 


DD(I1I1,I1,I.K-1) 
ЕЕ(11,11.1,К-1) 
ММ(11,11,1,К-1) 


DO 17 K 


GG(1,I,K-1 
66(2 , ] ,(-1 
66( 1 , [ ,(-1 
GG(4,1,K-1) 


CONT INUE 
CONTINUE 


1 , 4 

2 , KMAX - 1 

DTW / YACOB(I,K) * OELTAT(I.K) 
DD(I1,I1,I.,K-1) - OT1 г YACOB(I,K-1) 
ЕЕ(11,11,1,К-1) - DT1 © YACOB(I,K+1) 
1. + 2. * DTW « DELTAT(I,K) 

2 , KMAX - 1 

DQ1(I,K) 

DQ2(I,K) 

DQ3(I,K) 

004(1,К) 


PERFORM BLOCK TRIDIAGONAL MATRIX INVERSION FOR THE ENTIRE PLANE 


CALL MATRX2( IMAX,KMAX) 


DO 18 I = 2 , IMAX - 1 
DO 18 K = 2 , KM1 

001(1.к) = GG(1,I,K-1) 
DQ2(I.K) = GG(2,I,K-1) 
DQ3(1,K) = GG(3,1,K-1) 
DQ4(I,K) = GG(4,1,K-1) 


18 CONTINUE 


Css» UPDATE FLOW VARIABLES AT INTERIOR POINTS. 
967 CONTINUE 


RMAX = 0. 
RUMAX = 0. 
RVMAX ж Q. 
EMAX = 0. 
DO 995 K = 2 , KM1 
DO 19 I = 2 , IMI 
Q1(1,K) = Q1(I,K) + DQ1(I,K) * YACOB( I .K) 
Q2(I,K) = 02(1,К) + DQ2(I,K) » YACOB(I,K) 
Q3(I.K) = Q3(I,K) + 003(1,К) • YACOB(I.K) 
Q4(1,K) = 04(1,К) + 004(1,К) > YACOB(I,K) 
19 CONTINUE 
CI!!! !IDETC^MINE WHERE IN FLOW FIELO DENSITY IS CHANGING MOST RAPIDLY 
DO 935 I = 2 , [MAX - 1 
IF (RMAX.LT.ABS(DQ1(I,K)*eYACOB(I,K))) THEN 
IR = | 
KR = K 
END IF 
RMAX = AMAX1(RMAX,ABS(DQ1(I,K) * YACOB(I,K))) 
RUMAX = AMAX1(RUMAX,ABS(DQ2(I,K) = YACOB(I,K))) 
RVMAX = AMAX1(RVMAX,ABS(003(1.K) = YACOB(I,K))) 
EMAX = AMAX1(EMAX, ABS(DQ4( 1 .K) * YACOB( 1 ,K) ) ) 
995 CONTINUE 
С IF((ITN-1)/100*100.EQ. ( ITN-1)) WRITE (6.3002) 
С IF(ITN .EQ. 0) WRITE (6,3002) 
Cit! ISELECT INTERVAL AT WHICH OUTPUT OF RESIDUALS IS DESIRED 
C IF((ITN-1)/10»10.EQ. ( ITN-1)) WRITE (6,3001) RMAX,RUMAX ,RVMAX , 
С 1EMAX , IR, KR 
RETURN 


3002 FORMAT(//,4X, 'DRMAX',11X, 'DUMAX' , 11X, 'DVMAX' , 11X, 'DEMAX' , 10X, 
МА" ЗХ, КА") 
3001 ҒОКМАТ(4(Е14.8,2Х),215) 
ЕМО 
C 
Coss29253999229999292999999999999999929909$$229992992990999999999999992399$9$09992399929599952529259 
С 
SUBROUT INE MATRX1( IMAX , KMAX ) 
COMMON/TRID/DD(4,4,161,41) , MM(4,4,161,41), EE(4,4,161,41), 
166(4,161.41) 
СОММОМ/5СКАТ/А(4,4,161),НН(4,4,161,41),С(4,5,161) 
REAL MM 
REAL L11.L21,L31,L41,L22.L32,L42,L33,L43,L44 
2 JI DI .L3I (nel 


THIS SUBROUTINE PERFORMS THE BLOCK TRIDIAGONAL MATRIX IVERSION FOR 
AN ENTIRE PLANE DURING THE XI- SWEEP 


ОООО 


ОО 1 11 = 1, 4 

DO 1 K = 2 , KMAX - 1 

AI = 1. / MM(1,1,1,K) 

GG(I1.1.K) = СО(11,1,К) з Al 
HH(11,1,1,K) = EE(11,1,1,K) * Al 
HH(11,2,1,K) = EE(11,2,1,K) е АЈ 
HH(11,3,1,K) = EE(11,3,1,K) • АЛ 
HH(I1,4,1,K) к ЕЕ(11,4,1,К) * AI 

1 CONTINUE 


00 1000 I = 2 , IMAX — 2 
DO 5 пет, 4 
00 2 К = 2 . KMAX — 1 


C(11,1.K) = GG(I1,I,K) — DD(11,1,1,K) * GG(1,I—-1,K) 
1 — DD(11.2.1,K) * GG(2,I-1,K) 
2 - DD(11,3,I,K) * GG(3,I-1,K) 
У - 00(11,4,1,К) * GG(4,1-1,K) 
2 CONTINUE 
DO 5 l2=1 , 4 
005 К = 2 , KMAX - 1 
А(11,12,К)- ММ(11,12,1,К) - 00(11,1,1,К) * HH(1,12,1-1,K) 
1 — DD(I1,2.1,K) » НН(2,12,1-1,К) 
2 - 00(11,3,1,К) » PEE 
2 - 00(11,4,1,К) » НН(4,12,1-1,К 


C(11,12+1,K)= EE(I1,12,1,K) 
5 CONTINUE 

DO 3 K = 2 , KMAX - 1 

(11 = A(1,1,K) 

Са. / L11 

012 = А(1,2,К) • (11 

013 = А(1,3,К) • 111 

014 = А(1,4,К) • (11 

(21 = A(2,1,K) 

L31 = A(3,1,K) 

L41 = A(4,1,K) 

[22 = A(2.2.K) – 121 • 012 

L2I = 1. / L22 

023 = (A(2,3,K) — 121 * U13) г (21 

U24 = (A(2,4,K) — L21 * U14) • 121 

L32 = A(3,2,K) - L31 = U12 

L42 = A(4,2,K) - 141 е 012 

L33 = A(3,3,K) - L31 «• 013 — L32 « 023 

re. y s. 

034 = (А(3,4,К) — [31 • 014 - 132 е 124) # 131 

L43 = А(4,3,К) - [41 • 013 — [42 • 023 

[44 = А(4,4,К) — (41 • 014 — (42 • 024 - (43 ® U34 

(41 = 1. / (44 

C(1,1,K) = C(1,1,K) е» ЕП 

C(2,1,K) = (C(2,1,K) — 121 »* С(1,1,К)) « 12] 
C(3,1,K) = (C(3,1,K) — L31 * C(1.1.K) 
3 
. 


1 - 132 С(2,1,К)) » |151 
C(4 1,K) « (С(4,1,К) - 141 C(1,1,K) - 142 е С(2,1,К) 
1 – L43 • С(3,1,К)) • 141 
E К) = C(1,2,K) г 111 
C. К) = (C(2,2,K) — L21 © C(1,2,K)) < (21 
En „К) = (С(3,2 ,К) — Те ск) 
1 – [32 * C(2.2.K) ) з 131 
С(+.2,К) = (С(4,2,К) — 141 « С(1,2,К) - (42 е С(2,2,К) 
1 - [43 • С(3,2,К)) • [41 


C(1.3.K) = C(1.3.K) е 111 

С(2,3,К) = (С(2,3,К) – (21 в С.Ю) е ЕП 
GC3,3,.K) e (C(3.3,K) = t31- e C1 SK) 

1 - 132 « С(2,3.К)) е 131 


C(4,3,K) = (C(4,3,K) — (41 • С(1,3,К) - 142 з C(2,3,K) 
1 - 143 е C(3,3,K)) е 141 


С(1.4,к) « С(1.4,К) • (11 


Er = (C(2,4,K) - L21 • С(1,4,К)) • (21 

C(3,4,K) = (C(3,4,K) - L31 е C(1,4,K) 

1 = [32 • С(2,4,К)) • (31 

C(4,4,K) = (C(4,4,K) — 141 » С(1,4,К) - 142 г С(2,4,К) 

1 – [43 • C(3,4,K)) е 141 
C(1,5,K) = C(1,5,K) * L1I 

C(2,5,K) = (C(2.5,K) - 121 = C(1 5 kK)) s lia 

C(3,5,K) = (6(3.5.K) — £31 = GOs 

1 - (32 s C(2,5,K)) • (31 

C(4,5,K) = (C(4,5,K) — [41 • С(1,5,К) — (42 • C(2,5,K) 

1 – [43 • С(3,5,К)) • [4] 
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тк < K) — U34 


e C(4,1,K) 
СИ К) ж С(2,1,К) - 024 е C(4,1,K) 
1 - 023 * С(3‚,1,К 
C(1,1,K) = C(1,1,K) — 014 • Ir 
1 = И СЗТ К) = 012 • €(2,1,K) 
ек = C(3,2,K) — 034 » C(4,2,K) 
C(2,2,K) « С(2,2,К) - 024 « C(4,2,K) 
1 = 225% С(/3/2/К) 
C(1,2,K) = C(1,2,K) - U14 * MORET 
1 = ОВО 2 K) = U12 =< C(2,2 .K) 
Gs. 5.K) mC(3,3.K) — U34 е СО 3 К) 
C(2.3,K) = C(2,3,K) – 024 • С(4,3,К) 
1 - 023 » С(3,3,К) 
С(1,3.К)+ж С(1,3,К) = 014 = С, З.К) 
1 = 013 * С(3,3,К) - 012 а С(2,3,К) 
C(3,4,K) = RR - 034 = C(4,4,K) 
С(2,4,К) « С(2,4,К) - 024 • С(4,4,К) 
1 - 023 » С(3,4,К) 
C(1.4,K) к С(1,4,К) - 014 © C(4,4,K) 
1 = 015 е С(3,4,К) - 012 г С(2,4,К) 
C(3.5,K) = к - 024 « C(4,5,K) 
C(2,5,K) = C(2,5,K) – 024 èe С(4,5,К) 
1 - U23 e C(3,5,K) 
С(1,5,К) = С(1,5,К) — U14 = C(W,5,K) 
1 — U13 = C(3.5,K) — U12 е С(2.5 К) 
3 CONTINUE 
С 
00 6 11 = 1, 4 
00 9 K = 2 , КМАХ- 1 
9 GG(I1,I,K) = C(11,1,K) 
DO 6 I2 = 1 , 4 
DO 6 K = 2 , KMAX - 1 
HH(11,12,1,K) = C(11,12+1,K) 
6 CONTINUE 
10900 CONTINUE 
С 
С BACKWARD SUBSTITUTION 
С 
ОО 7 1 = IMAX – 3, 1 , - 1 
DO 7 11 ແ] , 4 
DO 7 K = 2 , KMAX - 1 
GG(11,1,K) = 66(11,1,К) - НН(11,1,1,К) • 66(1,1+1,К) 
1 - НН(11,2,1,К) е GG(2,1+1,K) 
2 - HH(11,3,1,K) * GG(3,I+1,K) 
3 - HH(11,4,1,K) * GG(4,1+1,K) 
7 CONTINUE 
RETURN 
END 


C 


Сова вое осо аа о оао ооо оо оса оао ооо ооо аз аз ооо ооо оооозо 


C 
SUBROUT INE MATRX2( IMAX , KMAX ) 
COMMON/TRID/DD(4,4,161,41),MM(4,4,161,41),EE(4,4,161,41), 
54 161,41) 
COMMON/SCRAT/A(4,4,161),HH(4,4,161,41),C(4,5,161) 

REAL MM 
REAL L11,L21,L31,L41,L22,L32,L42,L33,L43, L44 
2,L11, L21,L31,L41 


C 
С THIS SUBROUTINE PERFORMS THE BLOCK TRIDIAGONAL MATRIX IVERSION FOR 
C AN ENTIRE J=CONSTANT PLANE DURING THE ZETA- SWEEP 
€ 

DO 1 11 = 1, 4 

001 1 = 2 , IMAX - 1 

A] = 1. / MM(1.1.1,1) 

GG(11.1,1) = GG(I1,1I,1) * Al 

HH(11,1,1,1) = EE(11,1,1,1) + Al 
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HH(1142, 1,1 ТЕСИН СЕА 
НН(11,3,1,1) = ЕЕ(11,3,1,1) » AI 
HH(11,4,1,1) = EE(11,4,1,1) * AI 
1 CONTINUE 
DO 1000 K = 2 „ KMAX - 2 
00 5 I= 1 , 4 
DO 2 ] = 2 , IMAX — 
C(I1,1,1) = GG(I1,I,K) DD(I1,1,1,K) 


1 


| ' | 一 


DD(I1,2,I,K) 


e COT KI) 
* GG(2,I,K-1) 
СОЕ КИ 
GG(4,1,K-1) 


e HH(1,12,1,K-1) 
e HH(2,12,1,K-1) 
HH(3,12,1,K-1) 
“ НН(4,12,1,К-1) 


I) - L42 = C(2,1,1) 


‚I)) е“ 141 


- [45 » С(3,2,1)) = bal 


- 145 » С(3,5,1)) е 141 


- 145 » С(3,4,1)) е 141 


2 DD(I1,3,I,K) + 
3 - 00(11,4,1,K) e 
2 CONTINUE 
DO 5 12 = 1 , 4 
DO 5 | = 2 , IMAX = 1 
A(11,12,1)= M(11,12,1,K) - DD(11,1,]1,K) 
1 - 00(11,2,1,K) 
2 - 00(11,3,1,K) = 
3 — DD(I1,4,I,K) 
C(11,12+1,1)= €&(11,12,1,K) 
5 CONTINUE 
DO 3 I = 2 , IMAX - 1 
L11 = A(1,1,1) 
О Т А ШИ 
012 = A(1,2,1) 
U13 = A(1,3,1) * 
014 = А(1, 4,1) 
L21 = A(2,1,1) 
(31 = А(3,1,1) 
L41 = A(4.1,I) 
L22 = А(2.2,1) - (21 • U12 
[21 = 1. / 122 
023 = (А(2,3.1) - 121 * U13) * L2I 
024 ж (А(2,4,1) - 121» 014) » (21 
152 « А(3,2,1) - 131 ® U12 
L42 = А(4,2,1) — [41 » U12 
L33 = A(3,3,1) -— 131 » 113 - 132 • U23 
US] "wat [33 
U34 < (A(3,4,1) — 031 • 014 — L32 » U24) » L3I 
L43 = A(4,3,1) — 141 s U13 — 142 • 023 
(|44 « А(4,4,1) - 141 « 014 - 142 ເ“ 024 - 143 + U34 
L4] = 1. / 144 
C(1.1.1) m,C(1,1,1) В 
Gro qs (C(2, 1,1) - 121» С(1,1, ШТІ + L2I 
C(3,1,1) = (C(3,1,1) ຈ 441; 
1 - .32.4462.4.1))"“%31 
C(4,1,1) = (С(4,1,1) — (41 • С(1,1, 
1 — 143 » СОЗ. 1 
С(1,2.1) C(1. 2.1) И 
C(2,2,1) = (0(2,2,1) = 2r = C( T2) =I 
С(3,2,І ж“ (С(3,2,1) - (ЭЕМССОТ 2 1) 
1 - 152 =< С(2.2.1)) * [ЁЗ] 
С(4,2,1) = (C(4,2,I) — (41 • С(1,2,1) — 142 • С(2,2,1) 
1 
С(1,. 3 Т) ест, Зуи 
C(2,3,1) z (C(2,3,I) - L21 ““““25 
C(3.3 71) m (665 За = D33179 C 1 531) 
1 = 132 =< C(2 5 5) s< ແນ 
C(4,3,1) = (C(4,3,1) -— L41 • С(1,3,1) — 142 е С(2,3.1) 
1 
2(1,4,]) = С(1,4,1) • (11 
С(2,4,1) « (С(/2: Т) - (21 » C(1,4,I)) » L2I 
C(3,4,1) = (C(3,4,1) - L31 г С(1,4,1) 
1 - [32 • С(2,4,1)) • 131 
С(4,4,1) = (С(4,4,1) — [41 • С(1,4,1) – 142 • С(2, 4,1) 
С. ЕО s hl 
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02.51) = sn ә ҚЕНЕЭМӘРТӨТФІ)) „ 121 
ГБ - (С(2,5,1) - ІЗІ « C(1,5, I) 
1 SLSR CRS Opa as) 
eS T ES (C(4,5,1) — L41 • С(1,5.1) — (42 * C(2,5,1) 
1 – [43 • С(3,5,1)) • [4] 
C(3,1,1) BE - 034 » С(4,1,1) 
С(2,1,1) « С(2,1,1) - 024 • С(4, 1.1) 
1 – 023 • С(3,1,1) 
С(1,1,1) = С(1,1,1) = 014: е"С(4,1,1) 
1 SACA 1) =- 012 ® C(2 1 Г) 
(C (3 2.1) Е - 034 » С(4,2,1) 
92 2 1) С(2,2,1) - Џ24 е С(4, 2,1) 
1 223996055. 2-19 
291. 2.1) С(1,2,1) - Џ14 • С(4,2,1) 
1 = UJ е еп = Ру - 1012 е С(2 2 1) 
605 3,1) С(3,3,1) - 0934» С(4,53,1) 
C(2,3,1) = C(2,3,1) - 124 в C(4.3.1) 
1 — U23 = C(3,3,1) 
C(1,3,1) = C(1,3,1) — U14 e C(4,3,1) 
1 = 013% С(33,1) = 12 w C(2.3.1) 
C(3,4,1) C(3,4,1) — U34 e C(4,4,1) 
С(2,4,1) = С(2,4,1) — 024 * С(4,4,1) 
1 zuo. COS 4,1) 
901 .4,1) С(1,4,1) — U14 e C(4,4,1) 
1 - U13 * С(3,4,1) — U12 e C(2,4,1) 
c 5 5. 1) = C(3,5,1) - 034 * C(4 S, 1) 
С(2,5,1) = С(2,5,1) - Џ24 • С(4,5,1) 
1 - 023 8 С(3,5,1) 
951.5. 1) C(1,5,I) - 014» С(4, 5,1) 
1 = U13 »* С(3,5,1) = 012 = C(2,5,1) 
3 CONTINUE 
C 
DO 6 11 = 1,4 
DO 9 1 = 2, be ー 1 
9 GG(I1,I,K) “ЖЕГІ1-1,1) 
DO 6 12 ml, a 
DO 6 I = 2 , [МАХ 一 1 
HH(11,12,1,K) = C(11,12+1,1) 
6 CONTINUE 
1000 CONTINUE 
e 
C BACKWARD SUBSTITUTION 
C 
DO 7 K = KMAX — 5, 1 ,- 1 
ро 7 11 = 1 , 4 
DO 7 1 = 2 , IMAX - 1 
GG(11,1,K) = GG(11,1,K) - НН(11,1,1,К) е» GG(1,],K+1) 
1 — HH(I1,2,I,K) * GG(2,1,K+1) 
2 — HH(I1,3,I,K) * GG(3,I,K+1) 
3 - HH(I1,4,I,K) * GG(4,I,K+1) 
7 CONTINUE 
RETURN 
END 
C 
Сосо ооо ооо ооо ооо ооо оо ооо соо фо оо ооо ооо о со оз ооо ооо ооо ооо о ооо ооо осооо 
6 
SUBROUTINE METRIC 
COMMON/F I X/OMEGA 
COMMON/DGR ID/DT , IMAX,KMAX, ITEL, ITEU 
COMMON/GRID1/X(161,41),2(161,41) 
COMMON/GRID/YACOB(161,41) 
COMMON/MTRIX/XIX(161,41),XIZ2(161.41),ZETAX(161,41), 
+26142( 16] , 41) , 
1XIT(161,41),ZETAT(161,41) 
С 
Ceee SUBROUTINE METRIC COMPUTES THE METRICS IN BOTH DIRECTIONS AND 
C THE UNSTEADY COEFFICIENTS ETAT, ETC. 


DO 1999 K = 1 , KMAX 
DO 1990 I = 1 , IMAX 
XTAU = OMEGA « Z(I,K) 
YTAU = OMEGA = (-X(I,K)) 


Cees PRESENT SET UP IS FOR FLOW PAST AN AIRFOIL. 

E 

CI! I CENTRAL DIFFERENCES AT INTERIOR POINTS, TWO-POINT ONE-3.DED 
С!!! ! ıDIFFERENCES DOWNSTREAM, THREE-POINT AT OTHER OUTER BOUNDARIES 


C 


IF(I.EQ. 1.0R. I. EQ. IMAX) GO ТО 19 
XXI = .5 » (X(I*1,K)-X(I-1,K)) 
ZXI = .5 ຈ (Z(I+1,K)-Z(I-1,K)) 
со TO не 


10 IF(I.EQ.IMAX) GO TO 16 


XXI = 1.0 * (X(2,K) - X(1,K)) 
ZXI = 1.0 = (Z(2,K) - Z(1,K)) 
GO TO 15 


16 XXI = 1.@ * (X(IMAX,K) — X(IMAX-1,K)) 


ZXI = 1.0 = (Z(IMAX,K) — Z(IMAX-1,K)) 


15 CONTINUE 


IF(K.EQ.1.OR.K.EQ.KMAX) GO TO 17 
XZET = .5 *e(X(I,K+1)-X(I,K-1)) 

ZZET = .5 «(Z(I,K*1)-2(I,K-1)) 

GO TO 29 

IF(K.EQ.KMAX) GO TO 18 

XZET = 2. = X(I,2)-1.5 * X(I,1) - .5 = X(I,3) 
22ЕТ «2. «7(1,2) - 1.5. 2(1,1) - .5« 7(1,3) 
GO TO 29 


18 XZET = 1.5 в Х(І,КМАХ)-2.е Х(І,КМАХ-1)%.5еХх(І,КМАХ-2) 


ZZET = 1.5 • 2(1,КМАХ)—2.• 2(1,КМАХ—1)+.5•2(1,КМАХ—2) 


20 CONTINUE 


YACOBI = ХХІ » 22ЕТ - Х2ЕТ s» ZXI 

YACOB(I,K) = 1. / YACOBI 

XIX(I,K) = ZZET « YACOB(I,K) 

XIZ(I,K) = -XZET в ҮАСОВ(1,К) 

XTAU = OMEGA • 2(1,К) 

YTAU = — OMEGA « X(I,K) 

XIT(I,K) = — XIX(I.K) * XTAU ~ XIZ(I,K) * YTAU 
ZETAX(I,K) = —ZXI » YACOB(I,K) 

ZETAZ(I,K) « ХХІ » ҮАСОВ(1,К) 

ZETAT(I,K) « - ?ЕТАХ(І,К) * XTAU - ZETAZ(1,K) * YTAU 


1000 CONTINUE 


RETURN 
END 


Сееееееееоесееееееоеееевееоеоееееееееееееееевееооеееоеевевоеоеоеееевевесоееееееоееееевеееееееее 


С 


OOOO 


SUBROUTINE DISSIP 
COMMON/FLOW/Q1(161,41),02(161,41),03(161,41),04(161,41) 
СОММОМ/РЕКТВ/001(161,41),002(161,41),003(161,41),004(161,41) 
COMMON/DGRID/DT , IMAX, KMAX, ITEL, ITEU 

COMMON/GRID/YACOB( 161,41) 

COMMON /DAMP /WW , WW T 

DIMENSION P(161),EPS(161),DIS1(161,4),DIS2(161,4) 


THIS SUBROUTINE ADDS THE FOURTH ORDER DISSIPATION TERMS TO THE 
RIGHT HAND SIDE 


[М1 = IMAX — 1 
КМ1 = КМАХ - 1 
IM2 = IMAX — 2 
КМ2 = KMAX — 2 


DO 18 K=2 , KM1 
COMPUTE SWTICHING FUNCTION BASED ON SECOND DERIVATIVE OF PRESSURE 


DO 1 I = 1 , IMAX 


АНА САТИ) _Г0О2( 1 Кожа но] К)ч•е2)/(2 »ол(1,К))) 
] =1 , IMAX 

IP2 = I + 2 

IF(I.EQ.IM1) IP2 = IMAX 

IM2 = I — 2 

IF(I.EQ.2) IM2 = 1 

IPI = I + 1 

IF(I.EQ.IMAX) IP1 = IMAX 

IM = I - 1 

IF(I.EQ.1) IM = 1 

IF(I.EQ.1) IM2 = 1 

IF(I.EQ.IMAX) IP2 = IMAX 

EPS( ! ) = ABS(P( IP1 )-2. *P( 1)+P( IM) ) /ABS(P( IP1)+2 .*P( 1 )+P( IM) ) 
VOL = 2. /(YACOB(I,K)+YACOB(IP1,K)) 

VOL = 1: 

0151(1,1) = (QI(IP1,K)-Q1(I,K))*vOL 

0151(1,2) = I VOU 

DIS1(1,3) = (Q3(IP1,K)—Q3(I,K))=*VOL 

HP1 = Q4(IP1,K)+P(IP1) 

HP = Q4(1,K)+P(1) 

HM1 = Q4(IM,K) + P(IM) 

HP2 = Q4(IP2,K) + P(IP2) 

HPM = Q4(IM,K)+P( IM) 

01$1(1,4) = (HP1-HP ) *VOL 

0152(1,1) = (O1(IP2.K)-3.*(O1(IP1.,K)-Q1( 1 .K) )-G1( 1M,K ) ) VOL 
0152(1,2) —- (Q2(IP2.K)-3. *(Q2(IP1,K)-O2(1,K))-Q2( IM,K) )* VOL 
0152(1,3) = (QOS3(IP2,K)-3.*(Q3(IP1,K)-Q3(I1,K))-Q3( IM, K))*VOL 
0152(1,4) е- (НР2-3. » (НР1-НР)-НРМ) + VOL 

CONTINUE 

DO 15 1 = 1, IM 

D2P = AMAX1(EPS( 1 ) , EPS( 1+1 ) ) 

222 60. • 02Р 

СТИ АМАХ 1 (0.0, ( ] .-622)) 

022 C22 е WW/YACOB(I,K) + DT 

С11 = С11 = WW/YACOB(I,K) • ОТ 


CI!1! 'SWITCH ON SECOND-ORDER AND SWITCH OFF FOURTH-ORDER DISSIPATION 


C!!!!!IN VICINITY OF SHOCKS 


C 


16 
10 


) + C22 • 0151 
) + C22 “ 0151 
) 
) 


0151(1,1) = С11 • 0152(1 (1,1) 
(1,2) 
10272 DIS) IS 
(1,4) 


0151(1,2) = С11 • 0152(1, 
ВЕ! Оран Сл! » 0152(1, 
DIS1(1,4) = C11 е DIS2(T, 
CONT INUE 

0016 І «2, ІМІ 
001(1.К) = 01$1(1,1) - eH RM 
OQ2(I,K) » DIS1(I,2) - OIS1(1-1,2 
003(1,К) = DIS1(I,3) - າສ 
004(1,К) = 01$1(1,4) - 0151(1-1,4 
CONTINUE 

CONTINUE 

K DIRECTION 


1 
2 
3 
4) + С22 • 0151 


00 30] >2 , IM1 
wT= @.5 = DT « W / ҮАСОВ(1.2) 

W3 = 0.5 + DT » W / YACOB( 1 ,KM1 ) 

001(1 ,2) =WT* (Q1(I,1) – 2. • 01(1,2) + 01(1,3)) 


1+D01(1,2) 


002(1.2) «ИТ» (Q2(1,1) — 2. » 02(1,2) + 02(1,3)) 


1+DQ2( 1 ,2) 


DQ3(1,2) =МТ» (03(1,1) - 2. е 03(1,2) + 03(1.3)) 
1+005( ] ,2) 

004(1,2) =WTe (Q4(1,1) - 2. * Q4(1,2) + Q4(1,3)) 
1+004(1,2) 

WT= W3 

DO1(I,KM1) =WTe (Q1(1,KM2) – 2. • Q1(1,KM1) + Q1(],KMAX)) 
1+001(1,КМ1) 

DO2(1,KM1) =WTe (02(1,KM2) – 2. • Q2(1,KM1) + Q2(1,KMAX)) 
1+DQ2( 1 ,KM1 ) 
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DQ3(I,KM1) =WTe (03(1,KM2) — 2. © Q3(1,KM1) + Q3(1,KMAX)) 
1+003(1, КМ1) 
004( 1 , “ຟ1 ) «ИТ» (04(І,км2) - 2. * Q4(],KM1) + Q4(],KMAX)) 
1+DO4( 1 ,KM1 ) 
DO 35 K = 3 , KM2 
ите — DT * WW / YACOB( 1 ,K) 
DQ1(I,K) =WTs (Q1(1,K+2) — 4. © Q1(1,K+1) + 6. е 01( 
11,K) - 4. e Q01(1,K-1) + Q1(1,K-2))+D01(1,K) 
DQ2(I,K) =WT=*(Q2(I,K+2) — 4. = Q2(1,K+1) + 6. © Q2( 
11,K) - 4. e 02(1,K-1) + Q2(1,K-2))+002(1,K) 
003(1,К) «МТ«(03(1,К%2) - 4. © Q3(1,K+1) + 6. * Q3( 
11,K) — 4. + 03(1,K-1) + 03(1,K-2))+D03(1,K) 
DO4(I,K) zWT*(Q4(I,K*2) — 4. © Q4(1,K+1) + 6. * Q4( 
11,К) — 4. е О4(1,К—1) + 04(1,К—2) )+004 (1 „К) 
35 CONTINUE 
3@ CONTINUE 


RETURN 
END 
@ 
Ce ee の oe の 中 み の も もみ の ゆみ のみ の あみ ゆか の ゆ の ゆ の ゆあ の の ゆい ゆ の の の ゆい ゆ の ゆ ゆい ゆ みゆ あの 中 あゆ の ゆ の ゆあ ゆあ の の の ゆ ゆ いみ の みゆ ゆあ ゆみ ふみの あや もの あや も の の ゆい ゆ の の ゆ の ゆあ の ゆい ゆ も の ゆあ の の の の ゆ ゆ あの の ゆ ゆ の ゆず 
C 
SUBROUTINE WALLBC 
COMMON /SURF /PSUR (161) 
COMMON/GRID1/X(161,41),Z2(161,41) 
COMMON/PAR/GAMMA , REYREF , ALFA, ALFAT, REDFRE,AMINF,ALFAI 
COMMON/DGRID/DT., IMAX,KMAX, ITEL, ITEU 
COMMON/GRID/YACOB( 161,41) 
COMMON/DAMP /WWt , WW] 
COMMON/MTRIX/XIX(161,41).X12(161,41),ZETAX(161,41), 
+ZETAZ(161,41), 
1XIT(161,41),ZETAT(161,41) 
COMMON / F LOW/01(161,41),02(161,41),03(161,41),04(161,41) 
DIMENSION C1(4) 
DIMENSION A(2,2),RHS(2) 
C!!!! APPLY BOUNDARY CONDITIONS ON THE CUT AND THE AIRFOIL SURFACE 
DO 9 I=ITEU, IMAX 
11 = IMAX + 1 =- I 
О1(1,1) = .5 • (01(1,2)+01(11,2)) 
02(1,1) =. 5•(02(1,2)+02(11,2)) 
Q3(1,1) = .5 * (Q3(1,2) + Q3(11,2)) 
04(1,1) = .5 • (04(1,2)+04(11,2)) 
Q1(11,1)=01(1,1) 
Q2(11,1)=02(1,1) 
Q3(11,1)=03(1,1) 
04(11,1)=04(1,1) 


9 CONTINUE 
DO 1 дай ТЕС 'ITEU 
K = 3 
C1(1) = XIT(I,K) 
C1(2) = XIx(I.K) 
СО ЕО ОЕК) 
UCON3 = (Q2(1,K)*C1(2)4+03(1,K)*C1(3)) 
ТЕКО) 
K = 2 


С1(1) = XIT(I.K) 
C1(2) » XIX(I.K) 
C1(3) = xIZ(1,K) 
UCON2 « (Q2(I,K)*C1(2)4Q3(I,K)*C1(3)) 


1/01(1,K) 
RHS(1) = 2. + UCON2 - UCON3 - XIT(I,1) 
С FOR VISCOUS FLOWS SET UCON TO ZERO ALSO 


IF(REYREF.GT.@.) RHS(1) = — XIT(I,1) 
A(1,1) = XIX(1,1) 

A(1,2) = XIZ(I,1) 

А(2,1) = ХЕТАХ(1,1) 

А(2,2) = ZETAZ(1,1) 


C 


— 


RHS(2) = — ZETAT(I,1) 
TEMP1 = A(1,1) 


TEMP2 = A(1,2) 
ТЕМРЗ = А(2,1) 
ТЕМР4 = А(2,2) 
DEN = 1. /(TEMP1 • ТЕМР4 — TEMP2 • ТЕМРЗ) 


A(1,1) * A(2.2) * DEN 

А(1,2) = — ТЕМР2 • РЕМ 

А(2,1) = — ТЕМРЗ • ОЕМ 

A(2,2) = TEMP1 е ОЕМ 

2161 1) = 2: • Q1(162) - 01(1,3) 

02(1,1) « 01(1,1)е(А(1,1)еЕН5(1)%А(1,2) «КН5(2)) 
03(1,1) « 01(1,1) *(A(2.1)*sRHS( 1)+A( 2.2) *RHS (2 ) ) 
CONT INUE 

DO 18 I=ITEL ,ITEU 

а 

Ш-СЗСГІ 27/01( 1,2 

P2=( GAMMA-1 . ) *(Q4( ! ,2)- の 9.5*Q1( ! ,2) *(U2*U2+W2 *W2 ) ) 
05 92(1.2)/01(1.3) 

Ww>=03(1.3)/01(1,3) 
P3=(GAMMA-1.)+*(Q4(1,3)-0.5*01(1,3)+(UJeUI+WJeW3)) 
Р1=(4. •Р2—РЗ)/3. 

РЗИК ( 1 )= (САММА•Р1—1. )/(.7• АМ]МЕ •е2) 

ນ ະ 62( 1 .1)/01(1 ,1) 

W1=03(1,1)/01(1,1) 

04( ! ,1)=P1/( GAMMAー-1 . )+9.5*Q1( 1 ,1)*(U1*U1+W1 sW1 ) 
RETURN 

END 


Coosoo......... 0.000000... の いみ の あ の ふゆの の ゆあ ふみの も の ふも のみ の の の の も ゆ の の の の の ゆい ゆ 中 ゆ ゆみ の ゆ 上 の の の の の ふゆ ゆい ゆ も の の の の も の も の の あふ の あの の の あふ や の ゆあ の 


C 


OOO 


10 


SUBROUTINE STRESS(ITN,DALFA) 
COMMON/FLOW/Q1(161,41),02(161,41),03(*61,41),04(161,41) 
COMMON/DGRID/DT , IMAX , KMAX .ITEL,ITEU 
COMMON/GRID1/X(161,41),2(161,41) 

COMMON /PAR/GAMMA ,REYREF ,ALFA,ALFA1,REDFRE, AMINF ,ALFAI 
COMMON/PERTR/0Q1(161,41),002(161,41),003(161,41),004(161,41) 
COMMON /MUTUR/CMU( 161, 41) 

DIMENSION AA(161,41) 
1 .RH1(161 ) .RH2( 161 ) .RH3( 161 ) , RH4( 161) 

COMMON/ LOG IC/RSTRT,PITCH,PLUNGE 

LOGICAL RSTRT,PITCH,PLUNGE 

ລອດ ເຂ 2] ປ) / QVCI,J) 

Ham --95(1.J) / OI(T,J) 

THIS SUBROUTINE ADDS VISCOUS TERMS TO THE RIGHT HAND SIDE 
GOGM = GAMMA / (GAMMA — 1.) 

IF(ITN.GT.10.0R. (RSTRT)) CALL EDDY(DALFA) 

COMPUTE U AND V 

КМАХМ1 = КМАХ - 1 

IMAXM1 = IMAX - 1 

PR = 1. 

00 10 К = 1 , KMAX 

DO 10 I = 1 , IMAX 

E = Q4(1,K) / Q1(1,K) — @.5 * (U(I,K)**»24V(I,K)»92) 
АА(Т,К) = СОСМ • Е 


COMPUTE TXX,TXY AND VISCOUS DISSIPATION AT I - 1 / 2 


00 30 ເ > 2 , KMAXM1 

DO 20 I = 2 , IMAX 

UXI = U(I,K) - U(I-1,K) 

УХ] = V(I,K) - V(I-1,K) 

АХ! = AA(I,K) — AA(I-1,K) 

UZET= .25e(U(1,K+1)-U(1,K-1)+U(1-1,K+1)-U(1-1,K-1)) 
VZET= .25s(V(1,K+1)-VW(I ,K-1)+VW(I-1,K+1)-V( I-1,K-1)) 
AZET= .25*(AA(1,K+1)-AA(I,K-1)+AA(1-1,K+1)-AA(I-1,K-1)) 
XXI = X(I,K) - X(I-1,K) 
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20 


30 


ZXI = 2(1,К) - 2(]-1,%) 

XZET= .25 • (Х(1,К+1)—Х(1 ,К—1)+Х(1–1,К+1)-Х(1–1,К-1)) 
ZZET= .25 » (Z(I,K*1)-Z(I.K-1)4Z(1-1,K*1)-Z(1-1,K-1)) 
ХАС = XXI • 22ЕТ — 2Х1 = XZET 

ХАС = 1. / МАС 

XIX = ZZET • УАС 

ZETAX= — ZXI * YAC 

XIZ = —XZET * YAC 

ZETAZ= XXI * YAC 

CNM = .5 e (CMU(I,K) + CMU(I-1,K)) 


UX = UXI * XIX + UZET * ZETAX 

VX = VXI = XIX + VZET e ZETAX 

AX = AXI ® XIX + AZET » ZETAX 

UZ = UXI * XIZ + UZET з ZETAZ 

VZ = VXI © XIZ + VZET = ZETAZ 

AZ = AXI +. XIZ + AZET • ZETAZ 

TXX = -(-4. * UX + 2. * VZ) е СМ / 3. 

TXY = СММ • (02 + VX) 

TYY = -CNM / 3. = (-4. e VZ + 2. * UX) 

R4 = ((U(I,K)+U(I—-1,K))e*TXX+(V(I,K)+V(I—1,K))=*=TXY)e@.5 


+ CNM / PR/(GAMMA — 1.) • AX 
sé = ((U(I,K)+U(I-1.K))eTXY+(V(I,K)+V(I-1,K))*TYY)+0.5 
+ CNM / PR / (САМА - 1.) © AZ 


DEBUG 

TURN OFF ENRGY DISSIPATION AND DIFFUSION 
R4 = @. 
S4 = @. 
RH1(I) = 6. 
RH2 ( 】 ) = (XIX •» ТХХ + Х12 • TXY) / ТАС 
RH3(I) = (XIX = TXY + XIZ • ТҮҮ) / ТАС 
RH4(I) = (XIX + R4 + XIZ » 54) / YAC 


DO 30 I = 2 , 1МАХМ! 


DQ1(I,K) = DQ1(I,K) + ЕН1(1+1) — RH1(I) 

DQ2(I,K) = DQ2(I,K) + RH2(I+1) - RH2(I) 

DQ3(I,K) = DQ3(I,K) + RH3(I+1) - RH3(1) 

DQ4(I,K) = DQ4(I,K) + RH4(I+1) — RH4(1) 

CONTINUE 

IN THE Z DIRECTION 

ОО 70 I - 2, ІМАХМ1 

DO 69 K = 2 , KMAX 

UXI = .25 e (U(I+1,K)-U(I-1,K)+U(1+1,K-1)-U(I-1,K-1)) 
VXI 2 .25 » (V(I*1,K)-V(I-1,K)*V(I41,K-1)-V(I-1,K-1)) 
AXI = .25 = (AA(I+1,K)-AA(I-1,K)+AA(I+1,K-1)-AA(I-1.K-1)) 
XXI = .25 ® (Х(1+1,К)-Х(1-1,К)+Х(1+1,К-1)-Х(1-1,К-1)) 
ZXI = .25 = (Z(I+1,K)-Z(I-1,K)+Z(1I+1,K-1)-Z(I-1,K-1)) 
UZET = U(I,K) - U(I,K-1) 

VZET = V(I.K) - V(I,K-1) 

AZET = AA(I,K) — AA(I,K-1) 

XZET = X(I,K) - X(I,K-1) 

ZZET = Z(1,K) - Z(I,K-1) 

УАС = XXI » 22ЕТ - ZXI » XZET 

ЖАС = 1. / ТАС 

XIX = ZZET = YAC 


ZETAX= — ZXI * YAC 
ХІ2 = -Х2ЕТ * YAC 
ZETAZ= ХХІ » ТАС 


CNM = .5 © (CMU(I,K) + CMU(I,K-1)) 
UX = UXI » XIX + UZET • ZETAX 

VX = VXI = XIX + VZET e ZETAX 

AX = AXI » XIX + AZET * ZETAX 

UZ = UXI » XIZ + UZET e ZETAZ 

VZ = VXI = XIZ + VZET * ZETAZ 

AZ = AX] = XIZ + AZET » ZETAZ 

TXX = -(-4. * UX + 2. » VZ) * CNM / 3. 

TXY = CNM * (UZ + VX) 

TYY = -CNM / 3. + (—4. s VZ + 2. » UX) 

R4 = ((U(1,K)+U( I,K-1 ) ) *TXX+(V( 】,K)+V( 』,K-1 ) ) *TXY ) *9 . 5 
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1 4 CNM / PR/(GAMMA — 1.) e AX 
S4 = ((U(I,K)+U(1,K-1))*TxY+(V(I,K)+VW(I,K-1)).TYY)*0.5 


1 + СММ / PR / (GAMMA - 1.) « А? 

R4 = Q. 

54 = Ө. 

RH1(K) = Ө. 

M - (?ЕТАХ = TXX + ZETAZ * TXY) / МАС 

RH3(K) = (ZETAX e TXY + ZETAZ • ТҮҮ) / ҮАС 
69 RH4(K) = (ZETAX • К4 + 2ЕТА2 • 54) / МАС 


ро 78 К = 2 , КМАХМ! 

DQ1(I,K) = DQ1(I,K) + RH1(K+1) — i 

DO2(1,K) = DO2(1,K) + RH2(K+1) — RH2(K 

DO3(1,K) = DO3(I,K) + າ — RH3(K) 

004(1.К) = 004(1,К) + RH4(K+1) - RH4(K) 
70 CONTINUE 


RETURN 

END 
E 
Сееееееееееееоеееосееееееегеееееееоеееееегееееееооееоееегеоегеееееееегеевеовсвоввеевевее 
С 

SUBROUTINE LOAD(CPS,CL.CD,CM,ALFAS) 

COMMON/GRID1/X(161,41),Y(161,41) 

COMMON/DGRID/DT , IMAX, KMAX, ITEL, | ТЕЦ 

DIMENSION CPS(161) 


C 
С THIS SUBROUTINE COMPUTES THE INVISCID CONTRIBUTIONS 
C TO LOADS ON THE AIRFOIL SURFACE 
C 
IMAXM1 = IMAX - 1 
CL = 0. 
CD ແ 0. 
СМ = Ө. 
DO 400 I = ITEL , ITEU — 1 
XL = .5 • (Х(1,1)+Х(1+1,1)) 
YL = .5 • (У(1,1)+У(1+1,1)) 
ОХ = Х(1+1,1) - Х(1,1) 
DY = Y(1+1,1) = Y(1,1) 
CPA = CPS(I+1) * .5 + CPS(I) е 
DCL = CPA « (-DX) 
DCD = CPA œ DY 
CL = CL + DCL 
CD = CD + DCD 
400 CM = CM + DCD ® YL — DCL * XL 
С 
DCL = CL « COS(ALFAS) - CD * SIN(ALFAS ) 
CD = CL * SIN(ALFAS) + CD * COS(ALFAS ) 
CL = DCL 
RETURN 
END 
€ 


Сееееегееееоосееегееееееееоееоеооеееоееоеоееоеееееоеоеегеоеееоеееееоееееееооеееееооеееевее 
(5 

SUBROUTINE WRAP(II,JJ,XSING,YSING,XP,YP,SO.AQ,YO) 

DIMENSION S0(161,4),Y0(41,4),A0(161,4),XP(1),YP(1) 


THIS SUBROUTINE UNWRAPS THE AIRFOIL AND STORES THE UNWRAPPED 
SURFACE IN ARRAYS A0 AND S0. IT ALSO DETERMINES THE STRETCHING 
IN THE ETA DIRECTION. 


ООООО 


IMID = (II + 1) / 2 
DY = .8 / (ЈЈ - 2) 
DO 1 J=2 , JJ 
Y = FLOAT(J-2) * DY 
1 Y9(J,1) = 1.25 е Ү/(1.-Ү»Ү) 
YO(1 , 1) = - YO(3,1) 
P] = 4, e ATAN ( 1.) 
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ANGL =w PI + PI 

U = XP(1) — XSING 

V = YP(1) — YSING 

Џ = 1. 

У = Ө. 

IIM1 = II - 1 

002 1 = 1 , II 

X11 = XP(I) - XSING 

Y11 = YP(I) — YSING 

ANGL = ANGL + ATAN2((UsY11-Vex11),(Usx11+VeY11)) 
R = SQRT(X11*e2 + Ү11ее2) 


U = х11 
V =Y11 
R = SQRT(R) 


АО(1,1) = R COS(.5 * ANGL) 
2 50(1,1) = В © SIN(.5 е ANGL) 


C!!!!!IF OUTPUT OF UNWRAPPED COORDINATES IS DESIRED 

С WRITE (6,1000) 

С WRITE (6,2000) (1,A0(1,1),S0(1,1),1 = 1 , 11) 
RETURN 


1000 FORMAT(1X, 'UNWRAPPED COORDINATES IN THE TRANSFORMED PLANE ') 
2000 FORMAT(IS , 2F20.8) 
END 
C 
Ces622596966096969660266600$609$994945969665$9090960$9059699$6$606066059960600949696$60606$9996099099969029 
C 
SUBROUTINE TABINT(XP , YP, XSING, YSING,N) 
DIMENSION XP(161),YP(161),S0(161),A0(161) 
U = XP(1) - XSING 
V = YP(1) - YSING 
U = 1. 
V = 0. 
АМСЕ = 8. e ATAN(1.) 
DO 1 I = 1,N 
X11 = XP(I) - XSING 
Y11 = YP(I) - YSING 
ANGL = ANGL + ATAN2((UeY11-VsX11) .(UeX11+V*Y11 ) ) 


R = SORT(X11e0e2 + Y11 ss 2) 
U = X11 

ແ. 

R = SQRT(R) 


АО(1) = К * COS(ANGL * .5) 
1 50(1) = R «€ SIN(ANGL • .5) 
Ох -(АӨ(М)-АӨ(1))/96. 
AoST = AO(1) 
DO 3 I = 1 , 97 
ХХ = FLOAT(I-1) е DX + 4051 
CALL TAINT(AO,SO,XX, YY,N, 3, NER. MON) 
XP( 】 ) = XX * XX — YY e YY + XSING 
3 YP(I) = 2. * XX e YY + YSING 
RETURN 
END 
E 
Сегегеоееогеееееееееееоеоееоеовоеоеооооогеееееееееееееееееоеееогеееееееееееееееееееевее 
E 
SUBROUTINE ТАІМТ(ХТАВ,ҒТАВ,Х,ҒХ,М,К,МЕК,МОМ) 
DIMENSION XTAB(1),FTAB(1),T(10),C( 10) 


NASA — AMES SUBROUTINE FOR POLYNOMIAL INTERPOLATION 
OF A TABULATED FUNCTION 


O O O O 


IF(N-K) 1,1, 2 
1 NER = 2 

RETURN 
2 IF(K-9) 3,3,1 
3 ТЕ(МОМ) 4,4,5 
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5 IF(MON-2) 6,7,4 
4 J = 0 
NM1 = N — 1 
DO 8 I = 1 , NMI 
IF(XTAB(I) - XTAB(I+1)) 9,11,10 
11 NER = 3 
RETURN 
9 J = J-1 
GO TO 8 
19 J = J+1 
8 CONTINUE 
MON = 1 
IF(J) 12 , 6 , 6 
12 MON = 2 
7 DO 153 Т.М 
ТЕ(Х — ХТАВ(1)) 14,14,13 
14 Ј = 1 
СО ТО 18 
13 CONTINUE 
GO TO 15 
бро 16 1 1 , N 
[Е (Х-ХТАВ(Т)) 16, 17,17 
шъ | 
со то 18 
16 CONTINUE 
15 J = N 
18 J = J — (K+1) / 2 
ЈЕ(Ј) 19,19,20 
19 J = 1 
20 ຟ ແ J + K 
IF(M — N) 21,21,22 
22 Ј = Ј - 1 
GO TO 20 
21 KP1 = K + 
JSAVE = J 
26 DO 23 L = |, KPI 
C(L) = X — XTAB(J) 
T(L) = FTAB(J) 
23 J = 4+1 
DO 24 J = 1,K 
I = J+1 
ШБ оосор SD) 
= [+ 
ІЕ(І-КР1) 25,25,24 
24 CONTINUE 
FX = T(KP1) 
NER = 1 
RETURN 
END 


1 


ເກ 


2 
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SUBROUTINE SING(N2,N,X,Z,XLE,YLE,TEA,TES,XSING,YSING,CHD) 


THIS SUBROUTINE COMPUTES SINGULAR POINT LOCATIONS. 


O O O O O O O 


DIMENSION X(2) , Z(2) 
NÚ = N2 

№1 = N2 + 1 

N2 — 1 

X(N1 ) 

Z(N1 ) 

X(N2 ) 

Z(N2 ) 

X(N3) 

Z(N3) 

Х2 әә 2 - Х1 во 2 


2% 
кә 
MM ll uw n Ww H dq 
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D2 
03 
04 
05 
06 
07 
08 
В 


22 ее —27] 249 2 

(2 - X1 

22 - 21 

X3 we 2 ー X1 во 2 

23 se 2 - 21 ຂຂ 2 

X3 — X1 

23 - 21 

(07 • ( 01 + 02) – 03•(05+06))/(2. • (07•04–—03•08)) 


IF(ABS(D3).LT.ABS(D7)) GO TO 18 
А = (01 + 02 – 2. • В • 04) / (2. e D3) 
60 10 20 
10 А = (05 + 06 -2. • В • 08) / ( 2. = 07) 
20 CONTINUE 


R = 
XLE 
YLE 
CHD 


S 


ОКТ( (Х2-А)ее 2 + (22-В)ее2) 
X (NU) 
Z(NU) 
X(1) — X(NU) 


А2 - (2(2)-2(1)) /(Х(2) - Х(1)) 
A1 = (Z(N)-Z(N-1 ) )/(X(N)-X(N-1 ) ) 
TES = .5 » (АТ + А2) 


TEA 
TEA 


XSING 
YSING 
RETURN 


END 
С 


А2 - А1 

TEA » 57.29578 
= (A+XLE) /2. 
= (B+YLE) / 2. 


Ce oe みみ の の ゅ の の の の の の の の よ の の 上 ゆ の ゆ の の の の の ひひ の ひみ の ひよ ゆ の みよ ゆ の の ゆ の の みよ ひひ の の の ひ ゆ の の ゆみ ひよ の の の の よ ひ の ゆ ひひ の の の の の ゆ ひ の ゆ よみ よ の みゆ の の ゆ の の の ゆ ゆ の みよ の の の お の まず 


C 


SUBROUTINE AIRFOL(II,JJ,IT,IE) 
COMMON/GRID1/X(161,41),2(161,41) 

СОММОМ/Ү5ҮМ/ 15ҮМ 

DIMENSION S0(161,4),A0(161,4),Y0(41,4), XP( 161) , YP( 161), 
1ແ( 161 ) ,ແ(161) ,80(49) 


DATA (B0(1).1=1,32)/1.,1.0414,1.0836,1.1270,1.1715,1.2175,1.2651, 
11.3145,1.3659,1.4199,1.4755,1.5349,1.5973,1.6636,1.7342,1.8099, 
21.8914,1.9799,2.0764,2.1829,2.3012,2.4341,2.5653,2.7597,2.9646, 
33.2106,3.5141,3.9019,4.4219,5.1687,6.3632,8.6809/ 


C 


e s + + + 


008 [ແ] , 32 
8 А0(1,1) = 80(1) 

READ (5,1) 
1 FORMAT(1X) 

READ (5,2) FNU,FNL,ZSYM 
2 FORMAT(3F19.9) 

ISYM = а 


ЈЕ( 
II 


JJ 


eo е ва ël 
E 


25 


YM.NE.0.) ISYM = 1 
157 
41 
31 
127 


е АТАМ(1.) 


FNL 
NU * NL - 1 


READ(5, 1) 


READ (5,333) (XP(I),YP(I),I = NL , N) 
333 FORMAT(2F19.9) 
9994 FORMAT(F29.8) 
| = М + 1 
IF(ZSYM .GT. 0.) GO TO 9995 
L = NL + 1 
РЕАО (5,1) 
READ (5,333) (XP(L-I),YP(L-I),I=1,NL) 
GO TO 9996 
9995 K1 = L 
DO 16 I NL, N 
K = Ki 一 1] 
ХР(К) = ХР(1) 
УР(К) = - YP(I) 
16 CONTINUE 
9996 SCALE = 1. /(XP(1)-XP(NL)) 
XX = XP(NL) 
ҮҮ = YP(NL) 
DO 9997 I = 1 , N 
XP(I) = XP(I) з SCALE 
9997 УР(1) = YP(I) * SCALE 
CALL SING(NU,N, XP, YP, XLE, ZLE, TEA, TES, XSING, YSING,CHD) 
CALL TABINT(XP, YP, XSING, YSING,N) 
NBOCY = IE + 1 - IT 
DO 6791 I = 1 , NBODY 


1 
E(IT+L) = 

6791 F(IT+L) = 
IEP1 = IE 
SLOPT = 「 

DO 438 1 

II = I +1 — IE 

E(I) = А0(11,1) 

DXI = 1. / 48. 

D 240 sw (E(CI) 27.25) 

F(1) = F(IE) + SLOPT * ALOG(D) / D 


il 
=i 


L= IIP1 – 1 
Е(1) = Е(1) 
438 F(L) = F(IT) + SLOPT * ALOG(D)/D 
C WRITE (6,439) 
459 FORMAT(2X,3H 1,19Х,1НХ,19Х,1НҮ) 
C МАЛЕ ел и ЕСЕ О) аот 11) 
CALL WRAP(II,JJ,XSING,YSING,E.F,S@,AQ, YQ) 
C!!!!!MAP GRID BACK TO PHYSICAL PLANE AND SHIFT TO QUARTER CHORD 


DO 10 J "2 , JJ 
DO 10 I= 1 , II 
X( 1 .J-1 ) = АД(1,1)»•2 — (50(1,1)+У0(Ј,1) )••2 


] — 0.25 
10 2(1,)-1) = 2. * AQ(1,1) » (Se(1,1)+Y0(4,1)) 
е шу 
RETURN 
57 FORMAT( 15 .2F29.8) 
END 


С 
nenn 220 EEE EEE ERBE овоа оа ав ао сво を ゆり ваа оо ево 
し 
SUBROUTINE CLUSTR(DS) 
COMMON/GRID1/X(161,41),2(161,41) 
COMMON/DGRID/DT , IMAX,KMAX,ITEL,ITEU 
DIMENSION S(41),XP(41), YP(41),R(41) 
€ 


C THIS SUBROUTINE CLUSTERS A GIVEN X,Z GRID SUCH THAT THE FIRST POINT IS AT 
C THE USER-SPECIFIED DISTANCE DNMIN 
CI!!! COMPUTE THE OLD STRETCHING 

DO 100 I = 1 , IMAX 

S(1) = @. 


XP(1) = X(I,1) 
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YP(1) = Z(I ,1) 
DO 10 K = 2 , KMAX 
i - Х(І,к) 
ҮР(К) = 2(1,К) 
10 5(К) = SQRT((XP(K)-XP(K-1 ) ) **2+(YP(K )-YP(K-1 ) ) * * 2 ) 
1+S( K-1 ) 
SUMDX = S(KMAX ) 
CALL STRTCH( SUMDX , DS , F 1 .KMAX , FACTOR ) 
C WRITE (6,200) I,FACTOR 
R(1) = 0. 
DR = DS 
DO 20 K = 2 , KMAX 
R(K) = R(K-1) + DR 
DR = DR = FACTOR 
28 CONTINUE 
RLAST = 1. / R(KMAX) 
DO 30K = 2 , KMAX 
Ri = R(K) * RLAST » S(KMAX) 
CALL TAINT(S,XP,R1,XP1,KMAX, 3, NER, MON) 
X(I.K) = XP1 
CALL TAINT(S,YP,R1,YP1,KMAX,3,NER,MON) 
Z(I,K) = YP1 
30 CONTINUE 
100 CONTINUE 
С WRITE (6,115) 
DO 110 I = 1 , IMAX 
Ох + Х(1,2) - Se 
DY = 2(1,2) 20 33 


ОМ = SQRT(DXsDX+DYsDY) 
С WRITE(6,120) I , ОХ, DY, DN 
110 CONTINUE 
RETURN 


115 FORMAT(5X,6HNORMAL, 1X,8HDISTANCE,3H AT,4H THE,5H WALL, / 
1, 5Н 1, 8X, 2HDX,8X,2HDY, 8X, 2HON, //) 

128 FORMAT(15.3F10.5) 

299 FORMAT(I5,F19.5) 


END 
С 
Coooo..............0..0.......9.90......0..................................................... 
C 

SUBROUTINE STRTCH(SUMDX,DX1,F1,N1,R) 
C 
C THIS SUBROUTINE DETERMINES A GEOMETRIC 
С PROGRESSION OF GRID SPACING BETWEEN 1 AND N1 SUCH THAT 
С SUMBDX) EQUALS SUMDX. THE RATIO BETWEEN SUCCESSIVE 
E SPACINGS IS R. 

М = М1 - 1 

R= 1.5 

E1 = 1.E-94 

E2 = 1.E-04 

DO 19 L = 1, 50 

= (R-1) * SUMDX - ОХ1•(К••м—1) 

FP = SUMDX — DX1 * FLOAT(N) * R өв (М1) 

3...“ 
C IF(1.E-92.LT.RITER.AND.RITER.LT.1.) RITER = 1. 
C IF(1..LT.RITER.AND.RITER.LT.100.) RITER=.01 

IF(ABS(R-RITER) .LT. RsE1) GO TO 1 

R = RITER 

19 CONTINUE 

R= 1.0001 

Dx1 = DZTOT/FLOAT(N1-1 ) 

RETURN 

1 R= RITER 

RETURN 

END 
C 


12 
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С 


ооо 


OOOO 


2] 


SUBROUTINE EDDY(DALFA) 
COMMON/FLOW/Q1(161,41),02(161,41),03(161,41),04(161,41) 
COMMON/MUTUR/CMU( 161,41) 

COMMON/SK INCF/CF( 161) 

COMMON/DGRID/DT , IMAX,KMAX, ITEL,ITEU 

COMMON /PAR/GAMMA , REYREF , ALFA, ALFA1, REDFRE,AMINF , ALFAI 
ЕЕ отк an 

DIMENSION TIN(41),TOUT(41),Y(41 


INITIALIZE VISCOSITY EVERYWHERE 
FACTI = DT = AMINF / REYREF 
СМОМАХ = 100. າ FACT1 / DT 
DO 1 K = 1 , KMAX 
DO 1 | = 1 , [МАХ 
CMU(1,K) = FACTI 
THIS SUBROUTINE COMPUTES THE EDDY VISCOSITY BASED ON THE 
BALDWIN-LOMAX TWO LAYER MODEL 


DO 100 1 2 2 , IMAX 一 1 
UDIF = 0. 

FMAX = @.1E-06 

YMAX = .1E-06 


ҒҮМАХ = Ө. 
ү(1) = 0. 
UWALL = 0. 


IF(I.GT.ITEU.OR.I.LE.ITEL)UWALL = SQRT(Q2(I,1)**2+Q3(I,1)*e2)/ 


СЕ, 1) 


COMPUTE TAU AT THE WALL 


ОЕТ = 1.*(02(1,2)/01(1,2) — Q2(1,1)/Q1(1,1)) 

МЕ ແ າ 031 .2)/01(1;2) < Q3(I, 1)/Q1(1, 1)) 

XXI = Х(1+1,1) - X(I-1,1) 

ZXI = Z(I+1,1) - Z(1-1,1) 

ХЕТ = 4. ເ )(] 2) - 3. «® X(1,1) - X(1.3) 

ТЕТ = 4. * Z(I,2) — 3. | 2(1.1) - 2(1,3) 

ХХІ = .5 * XXI 

Хх [= „5 e ZXI 

XE 5 © ЛЕГ 

ТЕТ = 5 e ZET 

Aee (XXI + ZET = ZXI + XET) 

OMEGA = (ЧЕТ • ХХІ – VET » ZXI ) * YAC 

TwALL = AMINF © OMEGA / REYREF 

CF(I) = 2. * TWALL / (АМІМҒ» 2) 

FACT = SORT(Q1(1,1) » ABS(TwALL ) ) sREYREF/( 26 . *AMINF ) 
DO 19 K = 2 , KMAX-1 

UXI = (02(1+1,K)/01(1+1,K) - 02(1-1,K)/01(I-1,K)) 
VXI = (03(141,К)/01(1%1,К)-03(1-1,К)/01(1-1,к)) 
UET = (02(1,K+1)/01(1,K+1)-02(1,K-1)/01(1,K-1)) 
VET = (Q3(1.K+1)/Q1(1,K+1)-Q3(1,K-1) /Q1(1I,K-1)) 
XXI = X(14+1,K) - X(1-1,K) 

ZXI = Z(14+1,K) - Z(I-1,K) 

XET = X(1,K+1) — X(I,K-1) 

ТЕТ = 2(1,К+1) - 2(1,К-1) 

YAC = 1. / (XXI * ZET — ZXI • ХЕТ) 

OMEGA = ABS(UETeXXI+VETeZXI-UXIeXET-VXIeZET) * YAC 


UDIF = SORT(O2(1,K)+*2+03(1,K)=*2)/01(1,K) — UWALL 
IF(ABS(UDIF).GT.UDIFMAX) UDIFMAX = ABS(UDIF) 

Y(K) = SQRT( (X(I,K)-X(I,K-1) )**2+(Z( 1,K)-Z( 1 ,K-1) )**2)+Y(K-1) 
F = Y(K) ・ OMEGA 

IF((Y(K)eFACT).GT.20.) GO TO 31 

IF(I.GT.ITEL.ANO. I. LT. ITEU) F = F » (1. - EXP(-Y(K)«FACT)) 
CONT INUE 


MODIFIED TURBULENCE MODEL APPLY FOR SPECIFIED RANGE OF ANGLES WHERE 
FY IS USED TO FIND THE SECOND PEAK VALUE OF F FUNCTION 
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IF(ALFA.LT.ALFAI.AND.DALFA.GE.@.) THEN 
FY = F e Y(K) 
IF(FY.GT.FYMAX) THEN 
FYMAX = FY 
FMAX = F 
YMAX = Y(K) 
END I F 
ENDIF 
IF(ALFA.GE.ALFAI.OR.DALFA.LT.0.) THEN 
IF(F.GT.FMAX) THEN 
FMAX = F 
YMAX = Y(K) 
ENDIF 
ENDIF 
FCT = Y(K) * FACT 
IF(FCT.GT.20.) FCT = 20. 
FCT = ABS(FCT) 
Е = .4 * Y(K) е (1. - ЕХР(-ЕСТ)) 
TIN(K) = Q1(I,K) * EL * EL * OMEGA 
ТІМ(К) = ABS( TIN(K) ) 
19 CONTINUE 
KSwTCH = の 
FWAKE = YMAX e FMAX 
Fi = 0.25 = YMAX | UDIF «е2 / FMAX 
IF(F1.LT.FWAKE) FWAKE = Е1 
00 20 К = 2 , KMAX — 1 
FKLEB = 0. 
IF(ABS(Y(K)/YMAX ) . LT . 1.E+94 ) THEN 
FKLEB = 1. / (1. + 5.5 ® (09.3 e Y(K)/YMAX) »» 6) 
END IF 
TOUT(K) = .0168 «© 1.6 * Q1(I,K) e FWAKE * FKLEB 
TOUT(K) = ABS(TOUT(K) ) 
IF (KSWTCH.NE.@) GO ТО 20 
IF(TIN(K).GT.TOUT(K)) KSWTCH = K — 1 
20 CONTINUE 
DO 30 K = 2 , KMAX - 1 
IF(K.LE.KSWTCH) THEN 
CMU(I,K) = DT * (AMINF/REYREF + ABS(TIN(K))) 
ELSE 
CMU(I,K) = ОТ * (AMINF / REYREF + ABS(TOUT(K) ) ) 
END IF 
3@ CONTINUE 
100 CONTINUE 
RETURN 
END 
С 
Сзазоз ооо ао оз ос ооо оо соо о ооо о о ооо ооо оо соо ооо соо ооо о ооо зоо о оо засос оз аюсосо 
Е 
SUBROUTINE RESI(OMEGA) 
COMMON/PERTR/0Q1(161,41),DQ2(161,41),DQ3(161,41),D04(161,41) 
COMMON/GRID1/X(161,41),2(161,41) 
COMMON/DGRID/DT , IMAX, KMAX, ITEL, ITEU 
COMMON/FLOW/Q1(161,41),02(161,41),03(161,41),04(161,41) 
COMMON/PAR /GAMMA , REYREF , ALFA, ALFAT, REDFRE, AMINF , ALFAI 
DIMENSION RHS(161,4) 
XTAU(I,K) = OMEGA e Z(1,K) 
YTAU(I,K) = — OMEGA * X(I,K) 
THIS SUBROUTINE COMPUTES THE RESIDUAL ON THE RIGHT HAND 
SIDE ARISING FROM THE EULER- PART OF THE GOVERNING EQUATIONS 


O O O O 


FLUX TERMS WITHIN THE XI- DERIVATIVE 

DO 100 K = 2 , KMAX 一 1 

DO 1@ 1 = 1 , IMAX 

UCON = (Q2(1,K)/01(1,K) ) e (2(1,K+1)-Z(1,K-1)) 
1 - (Q3(1,K)/01(1,K)) e (X(1,K+1)-X(],K-1)) 

ОСОМ = 0.25» ЫТ» (СОМ 
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ХІТ = — XTAU(I,K) *(2(1.к+1)-2(1,К-1)) 
1 + YTAU(I,K) * (X(I,K*1) - X(I,K-1)) 
Mite Xi) е DT е 0.25 
UCON = UCON + XIT 
RHS(I,1) = UCON e Q1(1,K) 
ке! / QI(I,.K) 
P = (GAWA-1.) * (04(1,К) - .5 • R e(Q2(I,K)ee24 
1 ОЕ ал 
RHS(1,2) = Q2(I,K) + UCON + P e DT e 9.25 e (Z(1,K+1) — 2(1,К—1)) 
RHS(1,3) = O3(I,K) * UCON - P = DT e 0,25 e (X(1,K+1)-X(1,K-1)) 
RHS(I,4) = UCON = (Q4(1,K)+P) – ХІТ * P 
10 CONTINUE 
DO 11 】 = 2 . IMAX - 1 


001(1,К) = DQ!(I,K) - RHS(I+1,1) + RHS(I-1,1) 
002(1,К) = 002(1.К) - RHS(I+1,2) + RHS(I-1,2) 
003(1,К) = 003(1,К) - RHS(I+1,3) + RHS(I-1,3) 
11 DO4(I,K) = DO4(1,K) — RHS(I+1,4) + RHS(I-1,4) 
100 CONTINUE 
С 
С FLUX TERMS WITHIN THE ETA- DERIVATIVE 
С 


DO 200 1 - 2, ІМАХ - 1 
DO 20 К = 1 , КМАХ 
VCON = (02(1,K)/01(1,K) ) e (Z(I-1,K)-Z(I+1.K)) 
1 +(03(1,K)/01(1,K) ) e (X(I+1,K)-X(I-1,K)) 
VCON = VCON e 0.25 e DT 
ETAT = -XTAU(I,K) e (Z(1-1,K)-2(1+1,K)) — YTAU(I,K)e 
1 (X(1+1,K)-X(I-1,K)) 
ETAT = ETAT e 06.25 е ОТ 
VCON = VCON + ETAT 
RHS(K,1) = VCON e Q1(1,K) 
ວ > (GAMMA-1.) * (Q4(I,K) - 0.5 е#(02(1,К)ее2+03(1.К)ее2) /01 (Т,к)) 
RHS(K,2) = VCON * Q2(1,.K) +P © DT © .25 © (2(1-1,K)—Z(1+1,K)) 
RHS(K.3) = VCON * Q3(1,K) +P © DT © .25 e (X(1+1,K) - X(I-1,K)) 
RHS(K,4) = VCON • (04(1,К)+Р) — ETAT • Р 
20 CONTINUE 
DO 21 K = 2 , КМАХ - 1 
DQ1(I,K) = DQ1(I,K) - RHS(K+1,1) + RHS(K-1,1) 
DO2(1,K) = DO2(1,K) - RHS(K+1,2) + RHS(K-1,2) 
DO3(1,K) = 003(1,К) —- RHS(K+1,3) + RHS(K-1,3) 
21 DO4(I.K) = DO4(I,K) - RHS(K+1,4) + RHS(K-1,4) 
200 CONTINUE 
RETURN 
END 
C 
Ce*s5252202529209509022599006005060560559096095009009609060900600586006060560060060060909006000900090020900490900009952609560058909€ 
C 
SUBROUTINE ROTGRID(X,Z, IMAX,KMAX,DALFA) 
e ROTATE GRID IN THE CLOCKWISE DIRECTION BY AN AMOUNT DALFA 
DIMENSION X(161,41),2(161,41) 
CA = COS(DALFA) 
SA = - SIN(DALFA) 
DO 2@ K = 1 , KMAX 
DO 20 I = 1 , IMAX 
X1 = X(I,K) 
41 = 7(1.К) 
X(I,K) = X1 e CA - 21 е SA 
2@ 2(1,К) = Z1 e CA + X1 e SA 
RETURN 
END 


ооо оо ооо ооо оо соо со оо ооо осо ооо со ооо ос оо соо с ово во ооо ооо ово ооо ооо ово во ово 
SUBROUTINE CPPLOT(11,12,FMACH,X,Y,CP) 


THIS SUBROUTINE PLOTS CP AT EQUAL INTERVALS IN THE MAPPED PLANE 


OOO ООО 


По 


COMMON/SK INCF/CF (161) 
DIMENSION KODE(4).LINE(90),X(161),Y(161),CP(161) 
DIMENSION CFX(3,49),CFY(3, 49), CPX(3, 49) , CPY(3,49) 
DATA KODE/1H ,1H+,1H],1He/ 


в WRITE (6, 2) 
. 2 FORMAT(S@H@PLOT OF CP AT EQUAL INTERVALS IN THE MAPPED PLANE/ 
, 1 тено Х/С „лен СРЈ ,1ӨН CFU „лен СЕ. ) 


СРД = (1. + .2 • ЕМАСН ••2) го 3,5 — 1. 
600 = СРД / ( .7 • ЕМАСН ••2) 
КО = 30. • СРД + 4.5 
IMIN = (12-11)/2 + 11 
ILOW = 2 е IMIN 
ICOUNT = Ө 
снрех (11) - X(IMIN) 
00 12 1 = 1 , 90 
12 LINE(I) = KODE(1) 
DO 34 I = IMIN , 12 
= 30. = (CP® - CP(])) + 4.5 
КТ е 30. е (CP8 —P(ILOW-I)) + 4.5 
IF(K.LT.1) K = 1 
IF(K1.LT.1) K1 = 1 
IF(K.GT.98) K = 98 
IF(K1 .GT. 98) K1 = 98 
LINE(K0) = KODE(3) 
LINE(K) = KODE(2) 
LINE(K1) = KODE(4) 
ХОС = (Х(1) - X(IMIN)) / CHO 
+ WRITE (6,618) XOC,CP(I),CF(I),CF(ILOW-I).LINE 
LINE(K1) = KODE(1) 
34 LINE(K) = KODE(1) 
Cees GENERATE PLOT3D CP AND CF PLOTTING FILES 
ОО 500 ІкімІМ, 12 
XOC = (Х(1) - ДІНІҢ) D 
ICOUNT = ICOUNT + 
CPY(1,ICOUNT) = бого 
CFY(1,ICOUNT) » 0.000000 


CPY(2,ICOUNT) = CP(ILOW — I) 
CFY(2,ICOUNT) = CF(ILOW-1) 
CPY(3,ICOUNT) = CP(I) 

CFY(3, ICOUNT) = CF(I) 
CPX( 1 , ICOUNT ) = ХОС 
CFX( 1 , ICOUNT ) = ХОС 
CPX( 2 , ICOUNT ) = ХОС 

СЕХ(2, ICOUNT ) = ХОС 
CPX( 3, ICOUNT ) = ХОС 
CFX( さ 3, ICOUNT ) = ХОС 


500 CONTINUE 

IDM = 3 

JOM = 2 — IMIN + 1 

WRITE(58) IDM,JDM 

WRITE(SO)((CPX(1,J3), I=1,IDM),J=1,JDM), 
+ ((CPY(I,J), I=1,1DM),J=1,J0M) 
WRITE(69) 10M, JDM 

WRITE(69)((CFX( 1 ,J) , I=1,1DM),J=1.J0M), 


= ((CFY(1,J), I=1,IDM),J=1,JDM) 
RETURN 
610 FORMAT(4F10.4,90A1) 
END 


/EOF 


76 


TITELE: 
NACA 0012 AIRFOIL 


IMAX: KMAX: DT: WW: ALFA: ALFA1: ALFA]: REDFRE: AMINF: 
161 41 ‚005 Si 15.00 10.02 19. ее 0.151 .283 
ITEL: ITEU: REYREF: DNMIN: TSTART: FORMAT: RSTRT: PITCH: PLUNGE: 
5) 127 3.45 .00005 -1.0 3.0 TRUE TRUE FALSE 
СТР: CPLT: NSTP: PSTP: | 
0 0 1000 50 
FNU: FNL: FSYM: 
55. 33. js 

X: Y: 
0. 0. 
0.0050 „01153 

28125 .81894 

. 0250 02615 

. 0500 . 03555 

. 0758 04200 

. 1000 . 04683 

. 1500 ‚05345 

. 2009 2057357 

.2500 . 905941 

. 2600 ‚05962 

. 2700 . 05978 

. 2800 .85992 

. 2900 . 05999 

. 5000 06000 

. 3190 . 05999 

. 3200 ‚05992 

. 5500 05989 

. 5400 ‚05965 

‚ 5500 . 05947 

. 4000 .05803 

. 4500 . 25581 

. 5000 . 05294 

. 5500 04952 

. 6009 . 04563 

. 6500 041357 

.7000 . 03664 

. 7500 ‚03161 

. 8000 702022 

.8500 ‚02055 

. 9000 .01448 

.9500 .90807 
1.0000 .00126 
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